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Abstract 

V 

Nlany  computi 


Many  computational  problems  are  inherently  geometrical  in  nature.  For  example, 
cluster  analysis  involves  construction  of  convex  hulls  of  sets  of  points,  LSI  artwork 
analysis  requires  a  test  for  intersection  of  sets  of  line  segments,  computer  graphics 
involves  hidden  line  elimination,  and  even  linear  programming  can  be  expressed  in 
terms  of  intersection  of  half-spaces.  As  larger  geometric  problems  are  solved  on 
the  computer,  the  need  grows  for  faster  algorithms  to  solve  them.  The  topic  of  this 
thesis  is  the  use  of  geometric  transforms  as  algorithmic  tools  for  constructing  fast 
geometric  algorithms.  We  describe  several  geometric  problems  whose  solutions 
illustrate  the  use  of  geometric  transforms.  These  include  fast  algorithms  for 
intersecting  half-spaces,  constructing  Voronoi  diagrams,  and  computing  the 
Euclidean  diameter  of  a  set  of  points.  For  each  of  the  major  transforms  we  Include  a 
set  of  heuristics  to  enable  the  reader  to  usa  geometric  transforms  to  solve  his  own 
problems. 
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1 .  Introduction 

Many  computational  problems  are  inherently  geometrical  in  nature.  For  example, 
cluster  analysis  involves  construction  of  convex  hulls  of  points,  LSI  artwork  analysis 
requires  a  test  for  intersection  of  sets  of  line  segments,  computer  graphics  involves 
hidden  line  elimination,  and  even  linear  programming  can  be  expressed  in  'terms  of 
intersection  of  half-spaces.  As  larger  geometric  problems  are  solved  on  the 
computer,  grows  for  faster  algorithms  to  solve  them.  To  obtain  fast  geometric 
algorithms  a  set  of  tools  and  techniques  has  been  developed  that  takes  advantage 
of  the  structure  provided  by  the  geometry.  This  discipline  is  known  as 
Computational  Geometry.  In  the  following  sections  we  first  survey  the  previous  work 
in  computational  geometry  and  then  outline  the  contributions  of  this  thesis. 

1.1.  History  of  Computational  Geometry 

Geometry  has  been  studied  for  thousands  of  years  but  only  recently  has  it  been 
recast  in  computational  form.  Shames  [91]  describes  the  history  of  geometry  from 
the  perspective  of  a  computer  scientist.  Here  we  will  consider  only  the  history  of 
computational  geometry.  There  are  several  problem  areas  to  which  much  research 
has  been  devoted  —  construction  of  convex  hulls,  intersection  problems, 
closest-point  problems,  and  geometric  searching  problems.  In  this  section  we 
summarize  the  major  results  in  each  of  these  areas  and  also  a  number  of  topics  that 
do  not  fit  into  these  categories. 

1.1.1.  Convox  Hulls 

The  convex  hull  of  a  set  of  points  is  a  fundamental  geometrical  structure  that 
arises  in  a  multitude  of  different  problems  in  the  literature  and  this  thesis. 
Mathematics  texts  define  the  convex  hull  of  a  set  of  points  S  as  the  smallest 
convex  set  that  contains  all  of  the  points  of  S.  This  definition  is  fine  for  proving 
theorems  but  it  does  not  help  us  design  a  fast  algorithm. 

In  1970  Chand  and  Kapur  [25]  produced  a  convex  hull  algorithm  for  N  points  in 
K-space.  They  applied  a  procedure  called  "giftwrapping"  to  obtain  a  good  (but  not 
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Figure  1-1:  Convex  hull  of  a  planar  set  of  points. 

optimal)  algorithm.  Graham  [48]  presented  an  0(N  log  N)  time  planar  algorithm  in 
1972.  This  is  optimal,  in  the  worst-case  sense,  because  an  algorithm  that 
constructs  a  convex  hull  can  be  used  to  sort  [91].  (If  we  require  only  the  vertices 
of  the  convex  hull,  then  we  can  no  longer  use  a  convex  hull  algorithm  to  sort. 
Nevertheless,  Yao  [103]  has  proven  that  n(N  log  N)  time  is  still  required  in  the  worst 
case  when  only  quadratic  functions  of  the  input  are  allowed.)  Jarvis 
[54]  subsequently  applied  giftwrapping  to  the  planar  problem  to  obtain  an  O(VN) 
time  algorithm  where  V  is  the  number  of  vertices  on  the  hull.  If  V  is  less  than 
O(log  N)  then  Jarvis'  algorithm  is  faster  than  Graham's.  (Preparata  [81]  later  refined 
Graham's  result  by  constructing  an  0(N  log  N)  time  real-time  planar  convex  hull 
algorithm.  Rather  than  operating  on  all  N  points  collectively,  this  algorithm  updates 
the  hull  in  0(log  N)  time  after  each  point  is  read.)  Preparata  and  Hong  [83]  then 
solved  the  three-dimensional  convex  hull  problem  in  0(N  log  N)  time.  In  four 
dimensions,  however,  there  is  an  Q(N2)  lower  bound  because  the  convex  hull  can 
have  0(N2)  edges  ([49],  p.193). 


The  lower  bounds  above  apply  only  to  the  worst-case.  If  the  expected  number  of 
points  on  the  convex  hull  is  sublinear,  then  faster  expected-time  algorithms  are 
possible.  Floyd  [40]  and  Eddy  [39]  independently  discovered  a  planar  convex  hull 
algorithm  with  O(N)  expected-time  when  the  N  points  are  drawn  from  a  uniform 
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distribution  over  a  convex  region.  Bentley  and  Shamos  [16]  improved  that  result  to 
include  any  distribution  for  which  the  expected  number  of  points  on  the  hull  is  O(NP) 
for  some  p  <  1.  Furthermore,  their  result  extends  to  0(N)  expected-time  in 
three-space  while  maintaining  an  Q(N  log  N)  worst-case  time.  For  K  dimensions  we 
may  still  construct  the  convex  hull  in  0(N)  expected-time  if  the  K  coordinates  are 
drawn  from  independent  distributions.  In  this  case  the  expected  values  for  the 
number  of  maxima  and  the  square  of  the  number  of  maxima  are  only  O0og**^N)  [10] 
and  O(log^^'^N)  [30],  respectively.  Even  though  the  worst-case  number  of 
edges,  faces,  etc.  of  a  K-dimensional  convex  hull  of  N  vertices  grow  exponentially 
with  K,  the  expected  size  of  the  convex  hull  is  still  only  a  power  of  log  N. 

1.1.2.  Intersection  Problems 

Intersection  problems  have  also'  received  a  great  deal  of  attention  in  recent 
years.  They  occur  in  a  variety  of  areas  including  computer  graphics,  architectural 
data  bases,  printed  circuit  design,  and  even  linear  programming  [94,12,2,65]. 
Shamos  and  Hoey  [90,  91]  constructed  several  fundamental  intersection  algorithma 
including  a  linear  time  intersection  of  two  convex  polygons,  which  is  applied 
recursively  in  their  0(N  log  N)  time  algorithm  for  intersecting  N  half-planes.  Hoey*s 
0(N  log  N)  time  algorithm  [94]  to  determine  if  any  two  of  N  line  segments  intersect 
has  been  extended  by  Bentley  and  Ottmann  [12]  to  report  all  K  intersecting  pairs  in 
0(N  log  N  ♦  K  log  N)  time.  (Brown  [24]  has  reduced  the  storage  requirement  of 
Bentley  and  Ottmann's  algorithm  from  0(N  ♦  K)  to  0(N).)  If  ail  segments  are  either 
horizontal  or  vertical,  then  the  number  of  intersections  K.  can  be  counted  in 
0(N  log  N)  time  and  reported  in  0(N  log  N  ♦  K)  time.  (See  [4,  19,  98]  for  rectangle 
intersection  problems  and  algorithms.)  Finally,  Zolnowsky  [106]  and  Preparata  and 
Muller  [84]  have  applied  geometric  transforms  to  produce  three  algorithms  for 
intersecting  N  three-dimensional  half-spaces  in  0(N  log  N)  time. 
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1.1.3.  Closest  Point  Problems 

Closest  point  problems  arise  in  cluster  analysis,  pattern  recognition,  and,  in 
particular,  construction  of  minimum  spanning  trees  [93,  3].  Dobkin  and  Upton 
[34,.  35]  used  an  interesting  duality  transform  to  prove  an  Q(N  log  N)  time  lower 
bound  for  the  element-uniqueness  problem1  under  a  model  of  computation  that 
allows  (<,=,>)  comparisons  of  linear  functions  of  the  input.^  This  proves  an  fl(N  log 
N)  time  lower  bound  for  the  problem  of  finding  the  two  closest  of  N  points.  Their 
model  of  computation,  however,  allows  only  comparisons  between  linear  functions  of 
the  input.  With  a  stronger  model  of  computation  algorithms  faster  than  0(N  log  N) 
time  are  possible.  Fortune  and  Hopcroft  [41]  showed  that  if  the  floor  function  Is 
allowed,  the  two  closest  points  can  be  found  in  0(N  log  log  N)  time  in  the  worst 
case.  Previously,  Rabin  [86]  and  Yuval  [104]  had  given  O(N)  expected -time 
algorihms  for  the  K-dimensional  closest  pair  problem. 

Meanwhile,  Bentley  (and  others)  worked  on  multi-dimensional  nearest  neighbor 
problems  [9,  11,  14,  17]  and  he  invented  a  data  structure  called  a  K-D  tree  to 
solve  them  efficiently.  Bentley's  thesis  [3]  employed  a  strategy  called 
multi-dimensional  divide-and-conquer  with  which  he  obtained  the  first 
sub-quadratic  algorithms  for  several  multi-dimensional  closest  point  problems.  His 
thesis  is  also  a  good  source  for  learning  about  algorithm  design  —  rather  than  simply 
presenting  the  finished  product  he  displays  the  algorithm  design  process  and  at  the 
conclusion  presents  a  list  of  heuristics  to  use  in  designing  algorithms. 

Shamos  and  Hoey  [89,  91 ,  93]  created  on  0(N  log  N)  time  divide-and-conquer 
algorithm  for  constructing  a  Voronoi  diagram  of  N  planar  points.  A  Voronoi  diagram 

^The  element  uniqueness  problem  is  to  determine  whether  all  N  elements  of  an  unordered  multiset  are  unique. 

‘If  the  allowed  (<,=,>)  comparisons  arc  restricted  to  be  between  the  N  elements  themselves  —  no  linear 
functions  of  the  input  --  then  an  fl(N  tog  N)  time  lower  bound  is  easy  to  prove.  This  is  because  no  ordering  of  the 
N  eloments  other  than  a  total  ordering  can  guarantee  that  no  two  elements  are  equal.  Construction  of  a  total 
ordering,  however,  requires  a  sort,  which  costs  fl(N  log  M)  time. 
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(to  be  described  in  detail  in  Section  4.1)  contains  all  of  the  necessary  proximity 
information  to  solve  efficiently  a  number  of  closest-point  problems  including 
construction  of  a  Euclidean  minimal  spanning  tree,  a  proper  straighHine  triangulation 
of  the  N  points,  and  the  nearest  neighbor  problem.  Bentley,  Weide,  and  Yao 
[18]  hove  extended  the  techniques  of  Weide's  thesis  [99]  to  obtain  a  linear 
expected-time  algorithm  for  constructing  a  planar  Voronoi  diagram.  The  only 
conditions  are  that  the  probability  density  of  the  underlying  distribution  be  bounded 
both  above  by  a  constant  and  (below)  away  from  zero  over  some  finite  region. 

Bentley  and  Friedman  [8]  describe  a  heuristic  solution  for  a  minimum  spanning 
tree  algorithm  in  multi-dimensional  Euclidean  space  and  Yao  [102]  has  constructed 
provably  subquadratic  worst-case  MST  (and  related)  algorithms  for  the  Li,  Lg,  and 
Loo  metrics.  (For  K  >  3  dimensions  his  algorithms  take  0(N2*aOO(iog  N)1"«^)  time 
where  <*(K)  *  2"^*^.  For  the  special  case  of  the  Euclidean  metric  in  three 
dimensions  this  is  improved  to  0((N  log  N)1,8)  time.)  In  general,  though,  construction 
of  a  minimum  spanning  tree  is  a  graph  problem.  Kruskal  [63]  and  Prim  [85]  give  the 
classical  0(EV)  time  algorithms  (improved  to  0(V2)  time  [32]).  Both  Yao  [101]  and 
Cheriton  and  Tarjen  [26]  present  0(E  log  log  V)  time  algorithms  for  general  graphs  of 
E  edges  and  V  vertices,  and  Cheriton  and  Tarjan  [26]  also  present  an  O(V)  time 
algorithm  for  planar  graphs. 

1.1.4.  Geometric  Searching  Problems 

As  we  saw  above,  many  closest-point  problems  have  associated  searching 
problems;  in  this  section  we  summarize  a  separate  class  of  searching  problems  that 
are  not  related  to  any  particular  closest-point  problem. 

Shomos  [89]  gives  an  algorithm  that  locates  a  point  in  a  straight-line  planar  graph 
of  N  vertices  in  0(!og  N)  time  and  Dobkin  and  Lipton  [33]  and  Dewdney 
[31]  extended  this  technique  to  K  dimensions.  Unfortunately,  the  storage  and 
preprocessing  time  required  by  these  algorithms  are  prohibitive  —  0(N2)  in  the 

|Q 

planar  case  and  0(N2  )  in  the  K-dimensional  case  [102].  Lee  and  Preparata 
[66]  improved  the  storage  and  preprocessing  time  at  the  expense  of  increased 
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searching  time,  giving  an  algorithm  with  O(log^N)  query  time,  0(N  log  N) 
preprocessing  time,  and  0(N)  storage.  (Shamos  and  Hoey  [89,  93]  had  achieved 
these  bounds  for  searching  Voronoi  diagrams.)  Another  alternative  is  Preparata's 
algorithm  with  only  0(log  N)  query  time  but  0(N  log  N)  preprocessing  time  and  storage 
[80].  In  1977  Upton  and  Tarjan  achieved  an  0(log  N)  query  time  with  only  0(N  log  N) 
preprocessing  time  and  0(N)  storage  [69,  70].  The  point  location  problem  is 
generalized  to  the  location  of  a  sol  of  points  in  [68,  32]. 

Kung,  Luccio,  and  Preparata  [G4]  worked  on  the  problem  of  finding  the  maxima  of 
a  set  of  N  vectors  in  K-space.  (A  vector  is  maximal  if  none  of  the  other  N  -  1 
vectors  are  greater  in  all  K  coordinates.)  Using  divide-and-conquer,  they 
constructed  an  algorithm  that  finds  the  maxima  in  0(N  log  N)  time  in  two  dimensions 
and  0(N  (log  N)K"^)  time  in  K  >  3  dimensions.  They  also  proved  an  fl(N  log  N)  lower 
bound  for  the  problem.  Maxima  are  important  because  of  their  relationship  to 
convex  hulls  and  ECDF's  (to  be  described).  Bentley,  Kung,  Schkolnick,  and 
Thompson  [10]  extended  those  results  to  obtain  a  linear  expected  time  algorithm. 

Bentley  and  Shamos  have  created  a  fast  algorithm  for  constructing  and  searching 
an  empirical  cumulative  distribution  function  (ECDF)  [15].  An  ECDF  is  an  extension 
of  the  familiar  one-dimensional  cumulative  distribution  function.  The  value  of  the 
function  at  a  point  in  one  dimension  is  the  number  of  points  with  smaller  x 
coordinate.  In  K  dimensions  it  is  the  number  of  points  that  are  smaller  in  all  K 
coordinates.  Bentley  and  Shamos  applied  multi-dimensional  divide-and-conquer  to 
accomplish  ECDF  searching  in  0(logKN)  time  with  0(NlogK_1N)  storage  and 
0(N  logK-1N)  preprocessing  time.  Bentley  has  further  elaborated  this  in  his  papers 
on  range  searching  [6,  9,  11]. 

Bentley  and  Saxe  have  characterized  properties  of  a  large  class  of  problems 
called  decomposable  searching  problems  that  include  many  geometric  searching 
problems,  including  ECDF  searching  [5].  A  searching  problem  is  decomposable  if  the 
search  for  the  relation  of  an  object  x  to  a  set  S  =  A  U  B  satisfies 
query(x,  A  U  B)  =  query(x,  A)  *  query(x,  B) 
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for  any  sets  A  and  B  such  that  S  *  A  u  B  and  some  binary  function  M«M  that  is 
computable  in  0(1)  time.  Decomposable  problems  have  several  interesting 
properties,  one  of  which  is  that  any  static  searching  algorithm  for  a  decomposable 
problem  can  be  mechanically  transformed  to  a  dynamic  searching  algorithm  with  a 
loss  of  at  most  0(tog  N)  in  preprocessing  time  and  query  time.  They  have  extended 
this  result  to  a  class  of  alternate  preprocessing  time,  query  time,  and  storage 
tradeoffs  [13,  88]. 

1.1.5.  Miscellaneous  Geometric  Problems 

There  are  several  topics  that  do  not  properly  fit  into  any  of  the  above  categories. 
For  example,  we  should  mention  that  Garey,  Graham,  and  Johnson  and  Papadimitriou 
and  Steiglitz  were  the  first  to  demonstrate  that  several  geometric  problems  are 
NP-Complete  [45,  74].  Also,  Shamos  has  applied  many  of  the  techniques  of 
computational  geometry  to  statistics  and  created  a  new  field  of  computational 
statistics  [90].  ECOF  searching  (described  above)  grew  from  this  work,  and 
Weide's  thesis  [99]  gives  many  important  applications  of  statistical  techniques  to 
computer  science  problems.  This  includes  a  linear  expected-time  sorting  algorithm 
for  any  set  with  an  underlying  distribution  having  a  bounded  probability  density. 

1 .2.  Thesis  Outline 

The  topic  of  this  thesis  is  the  use  of  geometric  transforms  as  tools  for 
constructing  fast  geometric  algorithms.  The  object  of  using  a  transform  is,  of 
course,  to  give  the  problem  a  more  useful  representation  than  it  had  in  its  original 
form.  There  is,  however,  no  explicit  rule  for  determining  which,  if  any,  geometric 
transform(s)  can  be  applied  profitably  to  a  particular  problem.  Instead,  we  have 
generated  heuristics  for  application  of  the  geometric  transforms.  It  is  intended  that 
these  heuristics  will  help  the  reader  use  geometric  transforms  to  solve  his  own 
problems.  In  the  following  chapters  we  describe  a  collection  of  geometric  problems 
whose  solutions  illustrate  the  use  of  geometric  transforms.  The  algorithms  provide 
not  only  examples  of  the  applications  of  the  transforms  but  also  are  useful  results 
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by  themselves. 

The  use  of  transforms  is  not  new  to  computer  science.  For  example,  the  concept 
of  NP-completc  problems  (languages)  is  based  on  polynomial  time  reducibility  of  one 
problem  to  anotner  [1],  the  FFT  (Fast  Fourier  Transform)  is  used  for  fast 
multiplication  of  polynomials  [20],  and  many  "filter"  transforms  are  used  in  pattern 
recognition.  We  also  encounter  transforms  in  the  solution  of  difference  equations 
that  arise  in  the  analysis  of  algorithms  (the  2-transform  or  generating  function)  and 
solution  of  differential  equations  that  arise  in  analysis  of  several  types  of  queueing 
systems  (Laplace  transform)  [61].  Yet  another  common  example  is  obtaining  a  lower 
bound  on  the  complexity  of  a  problem  X  by  demonstrating  that  an  algorithm  that 
solves  X  can  solve  a  problem  Y  for  which  a  lower  bound  is  known.  (Shamos*  thesis 
[91]  gives  several  examples  of  this.)  Finally,  we  should  also  mention  Parker's 
thesis  [75],  which  explicitly  addresses  the  application  of  transforms  to  the  problems 
of  Huffman  tree  construction,  solution  of  nonlinear  recurrences,  and  construction  of 
permutation  networks. 

Chapter  2  establishes  the  context  and  direction  of  this  thesis  with  a  simple 
example  -  the  diameter  of  N  planar  points  in  the  L1  and  Lcq  metrics.  Before  this 
problem  is  discussed  we  first  settle  several  issues  --  model  of  computation, 
representation  of  the  geometrical  objects,  measures  of  complexity,  and,  of  course,  a 
definition  of  diameter  in  the  1-j  and  L&)  metrics.  We  present  an  algorithm  for 
computing  the  L,^  diameter  and  use  a  geometric  transform  (rotation)  to  transform 
the  L-j  diameter  problem  to  an  diameter  problem.  We  follow  the  same  schema 
followed  throughout  the  thesis;  Given  a  problem  Y  and  an  algorithm  that  solves  a 
(related)  problem  X,  we  apply  a  geometric  transform  f  that  transforms  problem  Y  to  a 
problem  of  type  X. 

In  Chapter  3  wc  describe  the  application  of  geometric  transforms  to  intersection 
and  union  problems.  We  solve  two  problems  in  detail  (the  intersection  of 
half-spaces  and  the  union  of  disks)  and  give  optimal  algorithms  for  each.  More 
importantly,  we  present  several  transforms  and  techniques  in  this  chapter  that  will 
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be  encountered  many  times  again  in  succeeding  chapters.  Hie  first  transform 
introduced  is  a  point  /  flat  duality  that  transforms  problems  that  involve  flats®  to 
(simpler)  problems  that  involve  points.  The  second  new  transform  (inversion) 
converts  problems  that  involve  circles  or  spheres  to  problems  of  lines  or  planes. 
Inversion  is  typically  combined  with  an  embedding  of  the  problem  in  a  higher 
dimension  to  add  another  degree  of  freedom  to  the  problem.  We  aiso  introduce  the 
convex  hull  (a  fundamental  geometrical  structure)  in  the  first  of  its  several 
applications  in  this  thesis.  Finally,  the  techniques  of  Chapter  3  are  tied  together  by 
deriving  the  point  /  flat  duality  from  a  limiting  case  of  inversion  combined  with  a 
convex  hull  and  a  linear  transform. 

In  Chapter  4  we  apply  inversion  (after  embedding  in  a  higher  dimension)  and 
convex  hulls  to  the  construction  of  nearest  point  tesselations  of  space.  The  most 
important  such  tesselation  is  the  Voronoi  diagram,  which  enables  efficient  solution  of 
a  number  of  geometric  problems  including  minimum  spanning  tree,  closest  points,  and 
Oelaunay  triangulation  of  a  set  of  points.  Shamos  [69,93,91]  applied 
divide-and-conquer  in  the  plane  to  obtain  the  first  0(N  log  N)  time  planar  Voronoi 
diagram  algorithm.  This  thesis  gives  a  new  0(N  tog  N)  time  algorithm  that  extends 
straighforwardly  to  higher  dimensions. 

Chapter  5  describes  two  surprising  applications  of  algorithms  that  search 
tesselations  and  the  transforms  used  are  the  same  in  both  cases  ~  the  point  /  flat 
duality  followed  by  an  orthographic  projection.  We  transform  linear  programming  in  K 
variables  and  N  constraints  to  a  problem  of  locating  a  point  in  a  K-1  -dimensional 
tesselation  induced  by  N  points.  The  problem  of  computing  the  Euclidean  diameter  of 
N  points  in  three-space  is  transformed  to  the  problem  of  finding  all  pairs  of 
overlapping  regions  in  two  outerplnnar  graphs  of  O(N)  vertices,  which  can  be  solved 
in  0((N  ♦  K)  log  N)  time  and  0(N)  storage  (where  K  is  the  number  of  pairs  of 
antipodal  vertices  of  the  convex  hut!  of  the  N  points). 


flat,  also  known  a*  a  hypcrplane,  prime,  or  a  (K-1)  flat,  it  a  K-l-ennentional  linearly  closed  tubspaee  of 
K-tpace.  Thus  a  line  it  a  flat  in  the  plane,  a  plane  it  a  flat  in  three-tpaee,  etc. 
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In  Chapter  6  we  cover  several  miscellaneous  problems  that  do  not  fall  in  any  of 
the  categories  of  the  previous  chapters.  We  describe  a  use  of  the  floor  function  to 
obtain  an  0(N  +  I/O  time  * -approximation  algorithm  for  the  Euclidean  diameter  of  N 
planar  points  and  also  demonstrate  an  application  of  gnomonic  projection  in  an 
algorithm  of  Yuval  [105]  for  determining  if  N  spherical  points  can  be  fit  in  a 
hemispherical  cap. 

Chapter  7  summarizes  the  thesis  and  points  out  directions  for  further  work. 
Appendix  I  describes  the  problem  of  choosing  a  good  orientation  for  flats  (before 
applying  the  point  /  flat  duality).  Appendix  il  gives  an  approach  toward  an 
0(N  log  N)  time  lower  bound  for  the  Euclidean  diameter  of  a  set  of  N  planar  points, 
and  Appendix  III  summarizes  the  geometric  transforms  used,  their  important 
properties,  and  their  applications. 
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2.  An  Example:  Diameter  in  the  Plane 

This  chapter  gives  some  of  the  flavor  of  geometric  transforms  by  presenting  a 
simple  example  ~  the  transformation  of  the  problem  of  computing  the  diameter  of 
a  set  of  N  planar  points  to  the  problem  of  computing  the  Lyj  diameter  of  another  set 
of  N  planar  points.  Out  first  we  require  some  definitions  and  explanations  along' the 
following  lines: 

1 .  a  precise  problem  specification, 

2.  a  model  of  computation, 

3.  an  appropriate  cost  measures  to  measure  the  complexity,  and 

4.  the  representation  of  the  problem  and  the  solution. 

In  the  following  sections  we  will  address  these  issues  and  then  construct  algorithms 
for  the  Li  and  diameters  of  a  set  of  planar  points. 

2.1.  Problem  Specification 

Let  S  *  {  P|  »  (Xj.yj),  i*1,...N  }  be  a  set  of  N  planar  points.  If  D(pf,pj)  equals  the 
distance  (in  some  as-yet  unspecified  metric)  between  points  Pj  and  pj,  then  the 
diameter  of  S  is 


DIAM(S)  =  D(Pj,p'). 

The  value  of  DIAM(S)  depends,  of  course,  on  the  metric  chosen  for  D.  The  three 
metrics  of  interest  in  most  applications  are 

h  metric:  D^pj.Pj)  *  |x,  -  Xj|  ♦  |yj  •  yj| 

L2  (Euclidean)  metric:  D2(Pj,Pj)  *  ( (xf  -  Xj)2  +  (yj  -  yj)2  )l/2 

L^  metric:  D00(pi,pj)  «  max(  |xj  -  XjJ,  (ys  -  yj|  ) 

Let  the  diameters  in  these  three  metrics  be  denoted  0IAM<](S),  DIAM2(S)  and 
OIAMco(S),  respectively.  The  unit  circles  for  these  metrics  are  pictured  In  Figure 
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Figure  2-1:  Unit  circles  in  the  L-j ,  I2,  and  L^o  metrics. 


2-1. 

The  problem  that  we  will  solve  is 

Given  an  algorithm  that  computes  DIAM^S)  for  any  set  S  of  N  planar 
points,  construct  an  algorithm  that  computes  DIAM^(S). 

The  solution  takes  advantage  of  a  natural  isometry  between  the  Lj  and  Lqq  metrics 

in  the  plane  [27].  We  will  cover  the  l_2  (Euclidean)  diameter  in  Section  5.2. 

2.2.  Mode!  of  Computation 

What  tools  are  we  given  to  compute  DIAMj(S)  and  DIAM^CS)?  In  other  words, 
what  is  the  model  of  computation?  There  are  two  (conflicting)  criteria  to  be  used  in 
our  choice:  (1)  how  realistically  the  model  reflects  the  capabilities  of  real 
machines,  and(2)  mathematical  tractibility  of  the  model.  The  real  RAM  [91]  (similar 
to  the  integer  RAM  [1])  is  a  reasonable  compromise  for  much  of  the  work  in 
geometric  algorithms.  Its  capabilities  are  basically  those  of  any  reasonable 
algebraic  programming  language  —  the  four  arithmetic  operations  (+,-,x,/), 
comparisons  between  numbers  (<,<),  and  indirect  addressing  (for  convenient  access 
to  arrays  and  other  structures).  A  word  in  a  real  RAM  is  assumed  to  be  able  to  store 
a  real  number  exactly;  although  this  assumption  is  not  entirely  realistic,  it  is  close 
enough  for  most  practical  applications.  We  often  augment  the  arithmetic  operations 
to  include  arbitrary  analytic  functions  (trigonometric  functions,  exponentials,  and 
logarithms,  etc.).  The  floor  function,  on  the  other  hand,  will  not  be  included  without 
special  comment  because  it  is  not  analytic. 
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The  floor  function  does  seem  to  add  power  to  our  model  of  computation  that  is  not 
available  from  analytic  functions  .alone.  Gonzalez  [47]  used  it  to  find  the  largest 
gap  between  N  (unsorted)  real  numbers  in  0(N)  time  and  Fortune  and  Hopcroft 
[41]  solved  the  closcst-point  problem  in  0(N  log  log  N)  worst-case  time.  Several 
fast  expected-time  algorithms  use  the  floor  function,  including  the  linear 
expectcd-timc  closest  point  algorithms  of  Rabin  [86]  and  Yuval  [104],  Weide 
[99]  uses  it  to  improve  his  linear  expected-time  sorting  algorithm  (for  all  underlying 
distributions  with  bounded  density)  and  Bentley,  Weide,  and  Yao  [18]  extend 
Weide's  result  to  linear  expected-time  Voronoi  diagrams  (for  certain  probability 
distributions). 

2.3.  Cost  Measures  and  Complexity 

Now  that  the  model  of  computation  has  been  defined  we  can  talk  about  the  cost 
or  complexity  of  an  algorithm  or  problem.  On  a  real  HAM  each  arithmetic  operation, 
comparison,  or  (indirect)  memory  reference  has  an  associated  cost.  The  cost  may  or 
may  not  depend  on  the  arguments  for  the  operation,  the  numbers  compared,  or  the 
contents  of  the  memory  referenced.  The  logarithmic  cost  criterion  for  an  integer 
BAM  [1]  does  assign  a  greater  cost  to  manipulations  (additions,  comparisons,  etc.) 
of  large  integers  than  for  small  integers.  But  for  a  real  RAM  it  makes  more  sense  to 
use  the  uniform  cost  criterion  —  all  operations,  comparisons,  and  memory 
references  have  a  unit  cost,  independent  of  the  numbers  being  manipulated.  We  will 
use  the  uniform  cost  criterion  throughout  the  thesis. 

The  cost  of  executing  an  algorithm  is  known  as  the  complexity  of  that  algorithm. 
The  complexity  of  a  problem  .is  the  minimum  complexity  of  any  possible  algorithm 
that  solves  it  (under  the  given  model  of  computation).  (The  complexity  of  an 
algorithm  is  always  an  upper  bound  for  the  complexity  of  the  problem  it  solves.)  The 
complexity  of  an  algorithm  or  a  problem  is  usually  expressed  as  a  function  of  the 
size  of  the  problem.  The  size  may  be  the  number  of  words  of  input,  output,  or 
whatever  is  most  appropriate  for  the  particular  problem.  It  is  often,  however, 
inconvenient  and  unnecessary  to  obtain  an  exact  count  of  ail  the  operations, 
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comparison,  and  memory  references  that  an  algorithm  makes  for  any  particular 
problem  size  N.  Knuth  [62]  has  popularized  a  convenient  notation  for  talking  about 
asymptotic  bounds  on  the  complexity  of  an  algorithm  or  problem: 

0(f(N))  =  set  of  all  functions  g(N)  such  that  for  some  positive 
constants  M  and  C,  |g(N)|  <  C  f(N),  for  all  N  >  M. 

fi(f(N))  =  set  of  all  functions  g(N)  such  that  for  some  positive 
constants  M  and  C,  g(N)  >  C  f(N),  for  all  N  >  M. 

0(f(N))  =  set  of  all  functions  g(N)  such  that  for  some  constants  M,  0-j, 
and  C2,  Ci  g(N)  <  f(N)  <  C2  g(N),  for  all  N  >  M. 

An  algorithm  that  solves  a  problem  of  size  N  in  f(M)  time  thus  proves  an  upper  bound 
of  0(f(N))  for  the  time  complexity  of  the  problem.  If  a  lower  bound  of  fl(f(N))  time  is 
also  known  for  that  problem,  then  that  problem  has  time  complexity  0(f(N)).  The 
complexity  of  an  algorithm  may  alternately  measure  the  space  or  storage  used.  The 
notation  is  the  same  as  for  time  complexity,  and  we  thus  may  speak  of  an  algorithm 
having  time  complexity  0(T(N))  and  space  complexity  0(S(N)). 

2.4.  Representation  of  the  Problem  and  Solution 

How  should  a  set  S  of  N  planar  points  be  represented  in  a  real  RAM  to  enable 
efficient  computation  of  DIAMOJ(S)  and  DIAMj(S)?  Many  data  structures  would  be 
suitable  but  the  simplest  is  either  an  N-by-2  array  or  two  arrays  X  and  Y  of  length  N. 
These  representations  are  reducible  to  each  other  in  linear  time.  Similarly,  different 
coordinate  systems  for  the  points  (X-Y  vs.  polar,  etc.)  are  linear-time  reducible. 
(The  solution  —  the  diameter  --  is  simply  a  scalar  real  so  its  representation  is  not  an 
important  issue  in  a  real  RAM.)  For  more  complicated  geometrical  objects  such  as 
polygons,  polyhedrons,  and  Voronoi  diagrams  the  issue  of  representation  is  not  as 
easily  solved,  and  those  problems  will  be  tackled  as  we  come  to  them. 
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2.5.  Algorithm  for  Lao  Diameter  of  a  Sot  of  Planar  Points 

The  Loj  diameter  of  a  set  of  planar  points  can  now  be  computed  fairly  easily. 
This  is  because,  as  shown  in  Figure  2-1,  the  circle  for  the  L^j  metric  is  a 
rectilinearly  oriented  square.  The  diameter  is  simply  the  diameter  of  the 
smallest  rectilinearly  oriented  square  that  contains  all  of  the  points.  The  diameter  is 
therefore  either  the  difference  in  y  coordinates  of  the  highest  and  lowest  points  or 
the  difference  in  x  coordinates  of  the  rightmost  and  leftmost  points.  Here  is  a 
pseudo-Algol  description  of  the  corresponding  algorithm: 

Lgo  Diameter  of  a  Set  of  Planar  Points 

Input:  integer  N  >  0,  arrays  X[1-.N]  and  Y[1:N] 

Output:  Lai  diameter  of  the  N  points 
Time:  O(N),  Storage:  Input  +  0(1) 

YMin  «■  YMax  +■  Y[1  ]; 

XMin  *•  XMax  X[1  ]; 
for  I 2  thru  N  do 
begin 

XMin  <-  min(  X[l],  XMin  ); 

XMax  ♦-  max(  X[l],  XMax  ); 

YMin  «-  min(  Y[l],  YMin  ); 

YMax  max(  Y[l],  YMax  ); 

end: 

L^Diameter «-  max(  XMax-XMin,  YMax-YMin  ) 

The  O(N)  time  complexity  of  the  above  algorithm  is  optimal  to  within  a  constant 
factor  because  the  algorithm  must  read  all  of  its  N  Inputs  to  ensure  a  correct 
answer.  There  is,  however,  room  for  improvement;  for  instance,  the  computation  of 
max  and  min  can  be  done  in  less  than  3/4  as  many  comparisons  as  are  taken  above 
[76].  Note  that  the  storage  required  is  actually  0(1)  rather  than  0(N)  because  no 
computation  involves  more  than  the  Ith  element  of  X  and  Y  at  any  given  time.  The 
values  in  X  and  Y  can  therefore  be  read  from  a  tape  rather  than  stored  in  arrays. 
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2.6.  Algorithm  for  L-j  Diameter  of  a  Set  of  Planar  Points 

In  this  section  we  will  construct  an  algorithm  for  the  L1  diameter  of  a  set  of  N 
points.  We  could  start  from  scratch,  but  since  the  La-,  diameter  algorithm  is  already 
available,  it  would  be  nice  to  be  able  to  make  the  L1  diameter  problem  look  like  an 
Loo  diameter  problem,  that  is,  transform  it  to  an  Lfx)  diameter  problem.  Fortunately, 
this  can  be  done,  and  the  clue  is  in  Figure  2-1.  The  circle  for  the  L1  metric  can  be 
made  to  look  like  the  circle  for  the  liX)  metric  if  it  is  simply  rotated  45  degrees  (and 
multiplied  by  a  scale  factor  of  2i/2.)  This  leads  us  to  the  intuitive  algorithm  below: 


L<]  Diameter  of  a  Set  of  Planar  Points 


Input:  integer  N  >  0,  arrays  X[1:N]  and  Y[1:N] 

Output:  L-j  Diameter 
Time:  O(N),  Storage:  O(N) 

!  Rotate  the  points  45  degrees; 
for  I  <-  1  thru  N  do 
begin 

X'  -  (X[l]  +  Y[l])  /  21/2; 

Y*  <-  (-X[l]  +  Y[l])  /  2i/2; 

X[l]«-X';  Y[l]  <-  Y'; 

end; 

!  Compute  l.u,Diameter  and  scale  by  2I/2; 

L1  Diameter  *-  DIAMt0(X,Y,N)  •  21/2; 

The  hard  port  is  proving  that  this  algorithm  is  correct.  Since  the  diameter  is 
simply  the  maximum  interpoint  distance,  it  will  be  sufficient  to  show  that  computing 
the  L-j  distance  by  the  definition  in  Section  2.1  is  equivalent  to  the  computation  in 
the  algorithm  above.  Let  p-t  =  (xj.yj)  and  Pj'  =  p;  rotated  n/4  radians  about  the  origin. 
(The  rotation  can  be  clockwise  or  counterclockwise,  as  long  as  it  is  the  same  in 
each  case.  In  the  algorithm  above,  we  use  the  formulas  x'  =  (x  +  y)  cos(lT/4)  and 
y'  *  (-x  +  y)  cos(n/4)  to  rotate  the  points  n/4  radians  clockwise.)  The  two 
methods  for  computing  the  L^  distance  between  Pj  and  pj  are: 
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1.  (Definition)  LjDistancetoj.Pj)  =  lx,  -  Xj|  +  |ys  -  yj|,  and 

2.  (Algorithm)  L1  Distance(Pj.Pj)  -  2l/z  *  LC0Distance(pj'>Pj'). 

We  prove  this  by  reducing  the  second  (algorithm)  formula  to  the  first  (definition) 
formula: 


V2  •  t(A)Oistancc(p1\pJ')  =  ft  •  max(  |x/  -  x/J.  jy,*  -  y/j  ) 


t-  r  x>y;  x,+y, 
=  J2  ■  max  -  -Lli 

l  72  J2 


5?ll 


33  3; 

V2  V2 

=  max(  ((x-x,)  +  (y-y,)!.  |-(xi-xj)  ♦  (yi-yJ)|  ) 
*  L,Distance(pi.pi) 


There  are  four  possible  cases  to  satisfy: 


1. x,  -  xj  <  o,  y|-yj<0, 

2.  x,  -  xj  2  0.  y,  -  yj  >  0. 


3.  x,  -  xj  1  0,  y,  -  yj  <  0,  and 


4.  x,  -  Xj  <  0,  yi-yj>0. 

in  each  of  these  cases  the  identity  holds.  The  algorithm  for  computing  the  L-, 
diameter  is  therefore  correct. 


2.7.  Principles  Covered 

In  this  chapter  we  solved  a  simple  geometric  problem  —  computing  the  Lf 
diameter  of  a  set  of  planar  points  —  and  demonstrated  the  use  of  a  geometric 
transform.  Several  principles  have  been  presented  that  will  be  encountered 
repeatedly  in  this  thesis:  precise  specification  of  the  problem,  choice  of  a  model  of 
computation,  cost  measures  and  analysis  of  the  complexity,  representation  of  the 
problem  and  its  solution,  and,  of  course,  the  use  of  a  geometric  transform.  The 
choice  of  the  transform  (rotation)  in  the  example  of  this  chapter  may  still  seem  like 
something  pulled  out  of  a  hat.  Yet  there  is  a  method  to  it,  as  the  following  chapters 
will  demonstrate. 


r 

i 
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3.  Intersection  and  Union  Problems 

In  this  chapter  we  introduce  two  important  techniques  —  the  use  of  a  point  /  flat 
duality  and  the  combined  use  of  inversion  with  embedding  in  a  higher  dimension  ~ 
and  apply  them  to  two  geometric  intersection  problems.  The  first  problem  is  the 
intersection  of  N  (UPPER)  half-spaces  and  the  second  is  the  union  of  (the  interiors 
of)  N  circles.  For  both  of  these  problems  we  develop  algorithms  that  are  optimal 
(within  a  constant  factor).  Finally,  the  last  section  of  this  chapter  shows  that  the 
techniques  that  we  used  for  these  two  problems  are  actually  more  closely  related 
titan  they  appear  to  be. 

3.1 .  Intersection  of  Half-Spaces 

In  this  section  we  analyze  the  problem  of  constructing  the  intersection  of  a  set  of 
N  (UPPER)  half-spaces.  The  first  topic  that  we  cover  is  the  representation  of 
half-spaces  and  their  common  intersection  in  a  computer.  Given  this  representation 
we  then  prove  upper  and  lower  bounds  on  the  complexity  of  constructing  the 
intersection  in  two,  three,  and  higher  dimensions.  We  conclude  with  fast 
expected-time  algorithms  and  some  open  problems. 

The  reader  should  carry  away  three  important  tools  for  the  construction  of 
geometric  algorithms: 

•  A  point  /  flat  duality  that  Is  applicable  to  a  number  of  problems  in  this 
thesis.  It  is  used  for  transforming  (formidable)  problems  that  involve 
flats  to  (simpler)  problems  that  involve  points. 

-  A  fast  algorithm  for  intersection  of  (UPPER)  half-spaces.  (An  algorithm 
for  intersection  of  half-spaces  is  useful  for  linear  programming  in  two 
or  three  variables  [91],  intersection  of  convex  polyhedra,  and  as  a 
tool  for  solving  other  geometric  problems.) 

•  The  first  of  several  important  uses  of  the  convex  hull  of  a  point  set. 

We  will  use  these  tools  many  times  in  succeeding  chapters. 
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3.1.1.  Representation  of  Half-Spaces  and  Their  Intersection 

The  first  requirement  of  any  representation  of  a  geometric  object  is  that  it 
contain  all  of  the  necessary  information  to  describe  the  object,  and  the  second 
requirement  is  that  it  provide  the  information  efficiently  (both  for  encoding  and 
decoding).  We  shall  first  describe  such  a  representation  for  the  two-dimensional 
case  (half-planes  and  intersections  of  half-planes)  and  then  extend  our 
representation  to  an  arbitrary  number  of  dimensions.  The  details  of  our 
representation  can  be  easily  modified  to  a  number  of  forms  that  can  be  reduced  to 
one  another  in  constant  time  for  a  single  object. 

We  represent  a  half-plane  by  the  line  bounding  the  half-plane  and  a  single  bit  to 
indicate  which  side  of  the  line  the  half-plane  is  on.  There  are  many  ways  to 
represent  the  boundary  line,  but  we  will  use  the  slope-intercept  form  with  the 
understanding  that  vertical  (or  near-vertical)  lines  will  require  exceptional  handling 
(Appendix  I).4  The  half-plane 

y  <  ax  +  b 

can  thus  be  represented  as  (a,b,0)  and  the  half-plane 

y  £  ax  +  b 

can  be  represented  as  (a.b.l).  In  a  computer  these  may  be  three  (scalar)  variables 
or,  if  there  are  many  half-planes,  three  elements  of  an  array. 

i 

The  representation  of  the  intersection  of  N  half-planes  is  more  interesting. 
Certainly  one  (cheap)  method  is  to  represent  the  N  half-planes  as  described  above 
and  include  a  scalar  flag  INT  that  indicates  that  the  intersection  is  intended.  This 
has  the  advantage  of  representing  the  intersection  fast  (in  linear  time)  but  the 
disadvantage  that  it  doesn't  help  us  quickly  answer  important  questions  about  the 
Intersection,  such  as  "Is  the  intersection  empty?".  Another  possibility  is  to  append 
to  each  of  the  N  half-planes  a  flag  that  indicates  whether  or  not  part  of  the 


^Preparata  and  Muller  [34]  use  a  homogeneous  coordinate  representation,  which  treats  all  coordinates  uniformly. 
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boundary  of  the  half-plane  is  also  a  boundary  of  the  intersection.  (If  the  boundary 
of  half-plane  I  doe3  not  meet  the  intersection  of  the  N  half-planes,  then  half-plahe  i 
is  redundant.)  This  representation  enables  us  to  answer  quickly  whether  or  not  the 
intersection  is  empty  but  in  the  worst  case  it  does  not  enable  a  faster  solution  to 
questions  of  the  form  "Is  point  P  inside  the  intersection?".  To  answer  such 
questions  quickly  we  must  store  the  nonredundant'  half-planes  in  sorted  order. 
Since  the  intersection  of  N  half-planes  is  a  (possibly  empty)  convex  polygonally 
bounded  region  with  at  most  N  edges,  the  representation  that  we  will  use  is  the 
quadruple 

(V,  M[1:V],  B[1:V],  F[1:V]). 

Here  V  is  the  number  of  edges  in  the  intersection,  M  and  B  are  the  slope  and 
intercept,  respectively,  of  the  lines  determined  by  the  V  edges,  sorted  in 
counterclockwise  order,  and  F  is  a  bit  vector  that  allows  us  to  quickly  distinguish 
the  inside  from  the  outside  of  the  intersection. 

The  representation  of  N  K-dimensional  half-spaces  is  a  simple  extension  of  the 
two-dimensional  case.  If  the  half-space  is 

k-1 

*  z 

i«1 

then  the  representation  is  simply 

(*ti  *2,  •  •  •  aK-1'  aK' 

Similarly,  if  the  "3"  is  replaced  by  a  ">"  in  the  equation  above,  then  the  "0"  wilt  be 
replaced  by  a  "1"  in  the  representation,  in  a  computer,  we  can  represent  N 
K-dimensional  half-spaces  in  one  large  array  A[1:N,0:K]  where  the  "side"  bits  are 
stored  in  the  entries  A[i,0]. 

%. 

The  intersection  of  N  K-dimensional  half-spaces  is  more  difficult  to  represent. 
This  is  because  the  total  number  of  vertices,  edges,  faces,  hyperfaces,  etc.  grows 
exponentially  with  K.  (If  we  choose  to  record  only  the  half-spaces  with  flags  that 
indicate  for  each  half-space  whether  or  not  it  is  redundant,  then  only  linear  storage 
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is  required.  Unfortunately,  as  we  mentioned  for  the  two-dimensional  case,  we  would 
not  then  be  able  to  answer  quickly  questions  of  the  form  “Is  point  P  in  the 
intersection?".)  We  must  first  establish  some  terminology.  Let  a  vertex  be  called  a 
0-face,  an  edge  a  1-face,  and,  in  general,  a  j-dimensionai  piece  of  the  intersection 
be  a  j-face.  We  will  represent  the  intersection  of  K-dimensional  faces  by 
enumerating  the  j-faces  (for  0  <  j  <  K-1)  and  recording  how  they  are 
interconnected.  The  representation  of  the  intersection  of  N  K-dimensional 
half-spaces  is  the  septuple 

(H,  F,  V,  Connect  1,  Connect  IPtr,  Connect2,  Connect2Ptr), 

where 

H[1:h,1:K]  represents  the  set  of  h  fiats  determined  by  the  K-1 -faces 
of  the  intersection, 

F[1:h]  is  a  bit-vector  of  flags  that  enable  us  to  distinguish  the  inside 
from  the  outside  of  the  intersection. 

V[1  :v,  1  :K]  is  the  set  of  v  vertices  determined  by  the  N  half-spaces, 

Connect  1  is  a  table  of  l-faces  used  by  ConnectlPtr, 

Connect1Ptr[l,J]  is  the  subscript  of  Connectl  for  the  first  1+1 -face 
that  the  Jth  l-face  bounds, 

Connect;?  is  a  table  of  l-faces  used  by  Connect2Ptr, 

Connect2Ptr[l,J]  is  the  subscript  of  Connect2  for  the  first  1-1 -face 
that  determine  the  Jth  l-face,  in  counterclockwise  order, 

Note  that  the  four  "Connect"  arrays  are  jagged  arrays  rather  than  rectangular 
arrays-.  They  are  also  redundant,  for  case  of  use. 


3.1.2.  Lower  Bound 

We  prove  an  ft(N  log  N)  time  lower  bound  for  the  intersection  of  N  haif-planes  by 
demonstrating  that  an  algorithm  that  intersects  half-planes  can  be  used  to  sort. 
(The  lower  bound  applies  for  all  half-spaces  of  dimension  K  £  2  because  half-planes 
are  just  a  special  case  of  K-dimensional  half-spaces.)  Our  construction  follows  that 
of  Shamos  [91]. 
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Theorem  1^  The  intersection  of  N  half-planes  requires  Q(N  log  N)  time  in 
the  worst  case. 

Proof:  Given  N  real  numbers  ajt  i=1,  ...  N  we  construct  N  half-planes  h| 
by 

hf:  y  £  ajX  *  (aj/2)2. 

These  hnlf-plnnes  hj  contain  the  origin  and  are  bounded  by  lines  that 
have  slope  aj  and  are  tangent  to  the  parabola  y  =  x2.  The  intersection  of 
the  hj  is  a  convex  polygonal  region  whose  edges  are  sorted  by  slope.  We 
simply  read  off  the  slopes  of  these  edges  to  obtain  the  aj  in  sorted  order. 

The  proof  of  the  lower  bound  for  intersection  of  half-planes  requires  a 
lower  bound  for  sorting  uiuier  a  model  of  computation  that  can  support  the 
operations  used  in  our  construction  above.  This  has  been  provided  by 
Friedman  [43],  who  proved  an  o(N  log  N)  time  lower  bound  for  sorting 
under  a  model  of  computation  that  allows  analytic  functions  of  the  input. 
Since  our  construction  requires  only  analytic  functions,  the  ft(N  log  N) 
lower  bound  for  sorting  applies  also  to  the  intersection  of  half-planes. 


3.1.3.  Intersection  of  Half-Planes 

Shamos  [89]  and  Shamos  and  Hoey  [94]  show  that  the  intersection  of  N 
half-planes  has  time  complexity  Q(N  log  N).  Their  algorithm  for  constructing  the 
intersection  in  0(N  log  N)  time  recursively  applies  their  linear-time  algorithm  for 
intersecting  two  convex  N-gons.  The  algorithm  that  we  describe  below,  on  the 
other  hand,  is  based  on  a  geometric  transform  (point  /  flat  duality)  that  maps  the 
problem  of  intersecting  half-planes  to  two  problems  of  constructing  the  convex  hull 
of  a  planar  set  of  points  (and  a  simple  intersection  problem.)  Furthermore,  it 
extends  to  higher  dimensions  (unlike  Shamos  and  Hoey's  algorithm).  We  next 
describe  the  decomposition  of  the  two-dimensional  problem  into  three  subproblems. 
The  following  sections  characterize  redundant  half-planes,  introduce  the  point  /  flat 
duality  transform  and  then  apply  the  transform  to  the  intersection  of  half-planes. 
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3.1. 3.1.  Tho  Two-Dimensional  Problem 

In  Figure  3-1  we  illustrate  the  intersection  of  N  half-planes.  The  intersection 
itself  is  indicated  by  the  shaded  region.  We  partition  the  half-planes  into  two  sets, 
UPPER  and  LOWER.  A  half-plane  is  in  set  UPPER  if  the  line  at  its  boundary  is  above 
the  rest  of  the  half-plane.  Similarly,  a  half-plane  is  in  the  set  LOWER  if  the  line  at 
its  boundary  is  below  the  rest  of  the  half-plane.  (If  any  boundary  lines  are  vertical, 
then  we  rotate  all  N  half-planes  a  small  angle.)  The  reason  that  we  produce  this 
partition  is  that  the  transform  (to  be  described)  actually  applies  only  to  lines,  not 
half-planes.  Since  each  line  may  be  associated  with  two  half-planes,  we  partition 
the  set  of  half-planes  into  two  parts  so  that  the  half-plane  determined  by  a  line  will 
not  be  ambiguous.  Our  partition  of  the  half-planes  also  enables  us  to  divide  the 
problem  of  intersecting  the  half-planes  into  three  parts: 

1 .  Construction  of  U,  the  intersection  of  the  UPPER  half-planes, 

2.  Construction  of  L.  the  intersection  of  the  LOWER  half-planes,  and 

3.  The  intersection  of  U  and  L. 

UPPER 


LOWER 

(a) 

Figure  3-1 :  (a)lnterseclion  of  N  half-planes.,  (b)  Intersection  of  regions  U  and  L. 

As  shown  in  Figure  3-1  b,  part  (3)  is  relatively  easy.  If  U  and  L  have  0(N) 
vertices,  then  the  intersection  (shaded  region)  can  be  constructed  in  0(N)  time.  We 


24  Oecember  1079. 


Geometric  Transforms 


PAGE  28 


describe  the  algorithm  in  detail  as  Algorithm  IntersectChains  below. 

Algorithm  IntersectChains 

Input:  Intersections  of  half-planes  U  *  (Ml,  UM[1:N1],  UB[1:N1])  and  L  *  (N2, 

LM[  1  :N2],  LB[1  :N2])  v/here  N1  and  M2  are  integers  such  that  N  *  N1  ♦  N2 

and  UM,  UB,  LM,  and  LB  are  the  slopes  and  intercepts  of  the  lines 

determined  by  the  edges  of  U  and  L.®  The  edges  are  sorted  in 

counterclockwise  order: 

UM[1  ]  <  UM[2]  <  ...  <  UM[N1] 

LM[1]  <  LM[2]  <  . . .  <  LM[N2] 

Output:  Integer  G  (number  of  vertices  of  the  intersection),  arrays  M[1:E],  B[1:E] 
(slopes  and  intercepts  of  the  E  edges),  bit  vector  F[1:E]  to  distinguish 
the  inside  from  the  outside  of  the  half-planes. 

Time:  O(N),  Space:  O(N). 


1 .  Scan  U  and  L  (vectors  UM,  UB,  LM,  and  LB)  from  left  to  right  until  two 
segments  intersect  at  a  point  P.  (If  no  segments  intersect  then  the 
intersection  of  U  and  L  is  empty.)  The  scan  can  be  done  in  O(N)  time  in 
a  manner  similar  to  the  0(N)  time  merge  in  the  merge  sort  algorithm. 

2.  Scan  U  and  L  from  right  to  left  until  two  segments  intersect  at  a  point 

Q. 

3.  If  P  ?  Q,  then  return  (in  vectors  M  and  B)  the  concatenation  of  the 
chains  of  line  segments  of  U  and  L  between  points  P  and  Q. 

4.  If  P  *  Q,  then  the  intersection  is  unbounded  (or  just  the  point  P  *  Q). 
In  the  case  of  an  unbounded  intersection  we  must  determine  whether 

return  the  chains  to  the  left  of  P  or  the  chains  to  the  right  of  P. 
Thjs  can  he  determined  by  comparing  the  slopes  of  the  rays  bounding 
the  left  and  right  sides  of  U  and  L.  If  the  slope  of  the  left  ray  of  U  is 
less  than  the  slope  of  the  left  ray  of  L,  then  return  the  chains  to  the 
left  of  P.  Otherwise,  return  the  chains  to  the  right  of  P. 


*Sinee  U  i»  an  interjection  of  UPPER  half -plane;  ana  L  is  an  interjection  of  COWER  half-plane*  it  i«  not 
necojjary  to  include  b<l  vector j  indicating  inside  vs.  outside  of  the  half -planes. 
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We  have  just  seen  iiow  to  construct  efficiently  the  intersection  of  N  half-planes, 
given  U  and  L,  the  intersections  of  the  UPPER  and  the  LOWER  half-planes.  Now  'we 
must  design  a  fast  algorithm  for  constructing  U  and  L.  Since  the  construction  of  L  is 
so  similar  to  the  construction  of  U,  we  wifi  describe  only  the  construction  of  U. 

3. 1.3. 2.  Redundant  Half-Planes 

Assume  that  the  N  half-planes  are  all  UPPER  half-planes.  Some  of  these 
half-planes,  such  as  half-plane  k  in  Figures  3-2  (a)  and  (b),  do  not  bound  any  side 
of  the  region  of  intersection.  It  would  be  nice  if  we  could  find  all  such  half-planes 
and  throw  them  away  since  they  do  not  contribute  to  the  final  result.  Once  that  is 
done,  we  can  find  the  intersection  of  the  UPPER  half-planes  rather  easily.  As  we 
can  see  in  Figure  3-1  (b),  the  slopes  of  the  sides  of  the  chain  U  are  monotonic 
decreasing  as  we  travel  from  left  to  right.  Given  the  lines  determined  by  the  sides, 
we  need  only  sort  the  lines  by  slope  to  determine  the  order  in  which  they  intersect 
to  form  the  sides.  Since  the  sort  costs  only  0(N  log  N)  time,  we  can  construct  the 
intersection  of  the  UPPER  half-planes  in  0(N  log  N)  time,  once  the  redundant 
half-planes  have  been  eliminated. 


Figure  3-2:  (a)  &  (b)  -  k  is  redundant,  (c)  -  k  is  nonredundant. 

How  do  we  determine  which  half-planes  are  redundant  and  which  are  not?  There 

are  two  conditions  that  we  need  to  check. 

Theorem  2a  An  UPPER  half-plane  k  is  redundant  with  respect  to  UPPER 
half-planes  i  and  j  iff 

(A)  Line  k  is  above  the  point  P  where  lines  i  and  j  meet,  and 
(D)  The  slope  of  line  k  lies  between  the  slopes  of  lines  i  and  j. 
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Proof:  The  half-spaces  i.  j,  and  k  are 

y  <  afx  +  bj,  y  <  ajx  +  bj,  and  y  £  akx  +  bk 

and  lines  i  and  j  meet  at  point 

P  =  (Px,Py)  =  ^  -fl>rbjMarap*  (ajbj-ajbjJ/Caj-aj)). 

(See  rirjnrr.  3-2.)  Half-plane  k  is  redundant  with  respect  to  half-planes  i 
and  ]  iff  line  k  lies  above  the  two  rays  rj  and  rj  originating  at  point  P  and 
defining  the  boundaries  of  the  intersection  of  half-spaces  i  and  j.  There 
are  six  possible  cases  to  consider: 

o  <  ai  <  aj-  0  <  aj  <  ah  ai  <  0  <  aj 

aj  <  0  <  ai-  ai  <  aj  <  and  aj  <  °i  <  °- 
Since  many  of  these  are  equivalent,  we  need  to  prove  only  the  two  cases 

0  <  aj  <  aj  and  aj  <  0  <  a,. 

Case  (1)  0  <  aj  <  aj:  The  ray  r{  points  (downward)  in  the  direction 
(1,-a,)  and  the  ray  rj  points  (upward)  in  the  direction  (1,aj).  Line  k  lies 
above  r(  and  rj  :ff 

Py  -  a;u  £  a^  (Px  -  u)  +  b^,  Yu>0  and  (1 ) 

Py  ♦  aju  S  ak  (Px  +  u)  *  b^.  Vu>0.  (2) 

Letting  u  *  0  in  either  (1)  or  (2),  we  see  that  line  k  lies  above  point  P, 
satisfying  condition  (A)  above.  To  prove  condition  (B),  that  lies 

between  aj  and  aj,  we  divide  by  u  in  (1)  and  (2)  and  then  take  the  limit 
asu-*  on,  obtaining  aj  >  ak  and  aj  <  a^,  respectively.  Conversely,  if  (A)  is 
satisfied  (Py  <  akPx  +  bk)  and  (8)  is  satisfied  (0  <  aj  <  a^  <  Bj),  then 
inequalities  (1)  and  (2)  immediately  follow. 

Case  (2)  aj  <  0  <  a^  The  ray  rs  points  (downward)  in  the  direction 
(1,-Oj)  and  the  ray  rj  points  (downward)  in  the  direction  (1,-aj).  Line  k 
lies  above  r(  and  rj  iff  inequalities  (1)  and  (2)  hold.  The  proof  is  very 
similar  to  the  proof  for  Case  (1).  □ 


How  fast  can  we  determine  (non)redundancy  for  each  of  the  N  UPPER 
half-planes?  Certainly  one  approach  is  to  test  all  pairs  of  half-planes  i  and  j  for 
each  half-plane  k.  That  costs  0(N®)  time,  though,  which  leaves  much  room  for 
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improvement,  in  the  next  section  we  show  an  entirely  different  way  of  looking  at 
this  problem  that  solves  it  in  a  natural  and  efficient  way. 


3. 1.3. 3.  A  Point  /  Line  Duality  Transform 

In  this  section  we  present  a  transform  that  exploits  a  natural  duality  between 
points  and  lines  in  the  plane.  A  line  in  slope-intercept  form  (y  =  ax  +  b)  is  uniquely 
identified  by  the  pair  (a,b).  (This  transform  will  not  work  for  vertical  lines.)  We  thus 
have  a  natural  mapping  from  liner,  to  points.  We  can  also  map  points  to  lines.  For  an 
arbitrary  point  (x,y)  the  set  of  all  lines  (in  slope-intercept  form)  that  pass  through 
that  point  can  be  represented  by  the  set  {  (a,b)  |  y  =  ax  +  b  }.  This  transforms  a 
point  (x,y)  to  a  line  b  =  -xa  +  y.  Points  thus  transform  to  lines  and  lines  transform  to 
points  by  the  formulas 

y  =  ax  +  b  -»  (a,b),  and 


(x.y)  ■*  b  =  -xa  +  y. 

This  duality  is  illustrated  in  Figure  3-3. 


Figure  3-3;  Point  /  Line  Duality 


This  transform  has  an  interesting  property:  Distances  in  the  y-coordinate 
between  points  and  lines  arc  preserved.®  The  difference  in  y  coordinate  between 


e 

The  restriction  in  tha  y  coordinate  is  important  because  it  can  be  snovn  to  be  impossible  to  preserve  the 
Suclldaan  distance  between  a  point  and  a  line  under  a  duality  transtotm  (02). 
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point  (c.d)  and  line  y  =  ex  *  f  is  d  -  (ec  ♦  f).  The  difference  in  the  transforms 
b  *  -ca  +  d  and  (e,f)  is  (-cc  +  d)  -  f;,  which  is  the  same.  It  follows  from  this  that 
incidence  is  preserved/  If  point  (c.d)  is  on  line  y  =  ex  ♦  f,  then  it  holds  also  for 
their  transforms  --  point  (c.f)  is  on  line  b  =  -ca  +  d.  Note  further  that  not  only  is  the 
magnitude  of  the  distance  (in  the  y-coordinate)  preserved  but  also  its  sign.  Thus, 
above/belownoss  is  preserved.  If  (c,d)  is  above  (below)  line  y  =  ex  ♦  f,  then  the 
transform  of  (c,d)  is  above  (below)  the  transform  of  y  =  ex  +  f. 

There  is  another  property  of  the  transform  that  we  should  mention.  The  transform 
is  not  involutory,  but  composition  of  it  four  times  produces  the  following: 

(x,y)  -»  b  =  -xa  +  y  -»  (-x,y)  -»  b  =  xa  ♦  y  -»  (x.y) 

Only  a  slight  change  is  required  to  make  the  transform  its  own  inverse:  express 
lines  in  the  form  y  +  ax  +  b  =  0  rather  than  y  =  ax  ♦  b.  Then  it  is  true  that 
y  +  ax  +  b  =  0«  (a,b).  But  this  Has  the  unfortunate  side  effect  that 
above/belowness  between  points  and  lines  is  not  preserved;  it  is  reversed.  If  point 
(c,d)  is  above  line  y  •*>  ex  +  f  =  0  then  the  transform  of  (c.d)  will  be  below  the 
transform  of  line  y  +  ex  +  f  =  0. 

3.1. 3.4.  Application  of  the  Transform  to  the  Two-dimensional  Problem 

We  now  show  how  the  transform  enables  us  to  intersect  the  UPPER  (or  LOWER) 
half-planes  fast.  More  specifically,  the  transform  enables  an  efficient  mechanism 
for  eliminating  the  "redundant"  half-planes.  Recall  the  two  conditions  for 
redundancy  of  an  UPPER  half-plane:  a  half-plane  k  is  redundant  iff  there  exist 
half-planes  i  and  j  such  that  (1)  line  k  is  above  the  point  P  where  lines  i  and  j 
intersect,  and  (2)  the  slope  of  k  is  between  the  slopes  of  lines  i  and  j.  In  the  ab 
plane  there  is  a  corresponding  interpretation. 

In  Figure  3-4,  line  k  is  above  point  P  in  the  xy  plane.  This  is  transformed  to  a 
point  k  that  is  above  line  P  in  the  ab  plane.  (Abovc/belowness  between  lines  and 


^ There  are  other  duality  transforms  that  preserve  incidence,  such  as  Ptuckcr's  transform  [91  ]. 
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Figure  3-4:  Transform  of  a  redundant  half-plane. 


points  is  preserved  by  the  transform.)  The  slope  of  a  line  in  the  xy  plane  is  the  a 
coordinate  of  the  corresponding  point  in  the  ab  plane.  Line  k  thus  has  a  slope 
between  the  slopes  of  lines  i  and  j  and  transforms  to  a  point  k  with  a  coordinate 
between  the  a  coordinates  of  points  i  and  j  in  the  ab  plane.  Figure  3-5  shows  the 
result  of  applying  the  transform  to  a  set  of  N  UPPER  half-planes.  A  point  in  the  ab 
plane  corresponds  to  a  redundant  half-plane  iff  it  is  directly  above  one  (or  more)  of 
the  line  segments  determined  by  the  N  -  1  other  points. 

— 


Figure  3-5:  Transform  of  N  UPPER  half-planes. 


Theorem  3:  Given  a  set  of  N  UPPER  half-planes  of  the  form 


and  a  mapping 


y  <  ajX  +  bj, 


y  <  ajX  +  bj  -»  (aj,  bj), 

the  nonredundant  half-planes  correspond  to  those  points  on  the  bottom 
part  of  the  convex  hull  of  the  N  points 


(af,  bj). 


Proof;  The  proof  is  in  two  parts:  (1)  a  point  on  the  bottom  part  of  the 
convex  hull  corresponds  to  a  nonredundant  half-plane,  and  (2)  a 
nonredundant  half-plane  transforms  to  a  point  on  the  bottom  part  of  the 
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convex  hull  of  the  N  points. 

1.  Let  P  be  a  point  on  the  bottom  part  of  the  convex  hull  of  the  N 
points  in  the  ab  plane.  P  does  not  lie  above  any  segment 
connecting  two  of  the  N  points  because  P  would  then  not  be  on 
the  bottom  part  of  the  hull.  It  follows  that  P  can  not  be 
redundant. 

2.  If  a  half-plane  P  is  nonredundant,  then  it  transforms  to  a  point  P 
that  docs  not  lie  above  any  segment  connecting  two  of  the  N  -  1 
other  points  in  the  ab  plane.  P  must  be  on  the  bottom  part  of  the 
convex  hull  because  otherwise  it  would  lie  above  such  a 
segment. 

□ 

We  have  reduced  the  problem  of  intersecting  N  UPPER  half-planes  to  the  problem 
of  constructing  the  (bottom  part  of  the)  convex  hull  of  N  points.  The  convex  hull  of 
N  points  in  the  plane  can  be  constructed  in  0(N  log  M)  time  [48].  This  leaves  only 
the  detail  of  separating  the  top  from  the  bottom  part  of  the  hull.  To  do  that,  we  find 
the  leftmost  point  of  the  hull  in  0(N)  time  and  then  traverse  the  hull  on  the  bottom 
side  until  the  rightmost  point  is  reached. 

Theorem  ^  The  intersection  of  N  half-planes  can  be  constructed  in 
0(N  log  N)  time. 

Proof:  We  have  broken  the  problem  of  intersecting  N  half-planes  into 
three  parts.  Part  (1),  the  intersection  of  the  UPPER  half-planes,  has 
been  shown  to  cost  only  0(N  log  N)  time.  Part  (2),  the  intersection  of  the 
LOWER  half-planes,  is  equivalent  to  part  (1)  so  it  can  also  be  done  in  0(N 
log  N)  lime.  Part  (3),  the  intersection  of  the  results  of  parts  <1)  and  (2), 
costs  only  O(N)  time.  The  intersection  of  N  half-planes  can  thus  be 
solved  in  0(N  log  N)  time.  □ 

It  is  interesting  to  note  that  we  can  also  use  an  intersection  of  half-planes 
algorithm  to  produce  a  convex  hull  of  points  algorithm.  We  first  transform  all  N 
points  to  UPPER  half-planes  by  the  formula 

(x,y)  -*  b  =  xa  ♦  y 
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and  intersect  the  half-planes.  We  transform  back  to  obtain  the  lower  part  of  the 
hull.  Then  we  transform  all  N  points  to  LOWER  half-planes  and  intersect  the 
half-planes  and  transform  back  for  the  upper  part  of  the  hull.  Merging  the  upper  and 
lower  parts  is  trivial,  since  the  leftmost  and  rightmost  points  will  be  in  each  one. 
The  total  time  required  is  0(N  log  N). 

3.1.4.  Intersection  of  Three-Dimensional  Half-Spaces 

The  technique  that  we  just  used  in  two  dimensions  can  be  extended  to  the 
intersection  of  N  UPPER  (or  N  LOWER)  three-dimensional  half-spaces.  (Zolnowsky 
[100]  and  Preporatn  and  Muller  [84]  describe  all  the  details  require  to  solve  the 
general  problem  of  intersecting  three-dimensional  half-spaces.)  The  following 
sections  extend  the  concept  of  redundant  half-space  to  three  dimensions  and  apply 
the  (three-dimensional)  point  /  flat  duality  to  construct  the  intersection  of  N  UPPER 
half -spaces  in  0(N  log  N)  time. 

2.1. 4.1.  Redundancy  in  Three  Dimensions 

The  algorithm  for  constructing  U  in  three  dimensions  is  analogous  to  the  algorithm 
for  the  two-dimensional  case.  We  first  transform  the  N  UPPER  half-spaces  to  N 
points  in  abc  space  and  construct  the  (bottom  part  of  the)  convex  hull  of  the  points 
in  0(N  log  N)  time.  The  points  above  the  bottom  part  of  the  convex  hull  correspond 
to  "redundant"  half-spaces  and  can  be  discarded.  To  form  U  wc  simply  apply  an 
inverse  transform  to  the  bottom  part  of  the  convex  hull.  We  will  now  describe  this 
procedure  in  detail. 

Assume  that  there  are  N  UPPER  half-spaces.  Some  of  these  half-spaces 
contribute  to  the  intersection  U  and  some  are  "redundant,"  such  as  half-space  M  In 
Figure  3-6.  For  the  plane  there  are  two  simple  conditions  for  redundancy  of  a 
half-plane.  In  three  dimensions  there  arc  two  analogous  conditions  for  a  half-space. 
The  first  condition  for  redundancy  of  a  half-space  M  is  that  plane  M  be  above  the 
point  P  where  planes  I,  J,  and  K  intersect,  as  in  figure  3-6.  The  second  condition, 
the  "betweenness  of  slopes"  condition,  is  more  complicated  to  express.  The 


it... . 
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purpose  of  the  "betweenness  of  slopes"  condition  is  to  ensure  that  a  plane  above 
the  point  P  cannot  drop  down  fast  enough  to  enter  the  region  below  planes  I,  J,  end 
K.  We  will  now  derive  an  algebraic  description  of  these  two  conditions. 

Let  plane  M  be  written 

2  3  «mx  ♦  »>My  ♦  CM 
and  let  planes  I,  J.  and  K  be  written 

2  =  a,x  ♦  b,y  ♦  c,, 

2  «  a0x  ♦  bjy  +  Cj,  and 
‘  2  *  aKx  ♦  bky  ♦  cK. 


The  point  P  *  (Px.  Py,  P2)  where  planes  I,  J,  and  K  meet  is  defined  by 

(3) 

The  three  rays  ry,  rjK,  and  r^j  can  be  expressed  as  vectors  originating  at  point  P. 
The  directions  of  these  three  vectors  are  obtained  by  computing  the  cross  products 
of  the  normals  to  planes  I,  J,  and  K.  For  example. 


-a,  -b, 

1 

rpx- 

'aj  *bj 

1 

py 

m 

cj 

.'•k  *bK 

1 

« 

px. 

«K. 
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I  J  K 

ai 

b|  -1 

ru  =  flu*  Yij)  = 

a,  b,  -1 

X 

aj 

bj  -1 

aj  bj  -1 

ak 

-1 

where  I,  J,  and  K  arc  unit  vectors  along  the  x,  y,  and  z  axes,  respectively. 
(Similarly  for  rj^  and  r«j.)  We  may  now  express  the  two  conditions  for  redundancy 
of  a  half-space. 

Theorem  5:  Let  UPPER  half-spaces  i,  J,  K,  and  M,  point  P,  and  rays  ry, 
rJK*  nncl  rKI  be  aG  Cl‘VRn  above.  Half-space  M  is  redundant  with  respect 
to  half-spaces  I,  J,  and  K  iff 

aMpx  +  bMpy  +  CM  -  pz 

and 

•m*m  +  bM<3,j  *  Yu 

a,V.«JK  +  Mjlv  *  7jk  (5) 

aM®Ki  +  brA.  *  Yki- 


Proof :  M  is  redundant  iff  it  lies  above  all  points  of  the  three  rays  ry, 
rjK,  and  rK).  But  M  lies  above  these  three  rays  iff 

<Px  +  U«IJ>  +  bM<Py  +  u/U  +  CN1  *  (pt  ♦  »?«)>  V">0. 

aM  <PX  +  U«JK)  +  bM  <Py  +  U0Jk>  +  CU  *  (Ps  +  «7jk).  V">n- 

aM(Px  +  u«u)  +  bM(Py  +  u|5H)  ♦  cM  >  (P4  +  "rw).  v">n 

Letting  u  =  0  we  obtain  the  first  condition  (inequality  (4))  for  redundancy 
of  M 

aMpx  +  bMpy  +  CM  -  Pz 

and  dividing  by  u  and  taking  the  limit  as  u  -*  ro  we  obtain  the  other  three 
conditions  (5).  Conversely,  if  the  four  inequalities  of  (4)  and  (5)  are 
satisfied,  then  M  must  lie  above  all  points  of  the  three  rays  ry,  r^,  and 

rKj.  We  prove  this  by  simply  multiplying  both  sides  of  the  last  three 
inequalities  (5)  by  a  u  >  0  and  adding  the  result  to  the  first  inequality  (4). 

□. 


m/mm  m  -  - 
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3. 1.4.2.  Application  of  the  Transform  to  ths  Intersection  of  Half-Spaces 

The  transform  that  we  use  in  three  dimensions  is  a  straightforward  extension  of 
the  two-dimensional  transform;  planes  transform  to  points  and  points  transform  to 
planes.  The  formulas  for  the  transform  are:® 

2  *  ax  +  by  +  c  -*  (a,b,c),  and  ^ 

(x,y,z)  -*  c  =  -xa  ♦  -yb  +  z. 


The  distance  between  a  point  and  a  plane  In  the  z  coordinate  is  preserved  by 

this  transform  and,  most  importantly,  the  sense  of  nbove/belowness  (in  the  z 

coordinate)  between  points  and  planes  is  also  preserved. 

Theorem  6:  Let  the  UPPCR  half-spaces  I,  J.  K,  and  M  transform  to 
points  P|,  Pj,  PK,  and  P^  by  the  transform  (G).  Half-space  M  is  redundant 
with  respect  to  half-spaces  I,  J,  and  K  iff  point  PM  is  directly  above  some 
point  in  (or  on)  the  triangle  determined  by  points  P|,  Pj,  and  PK. 

Proof:  Plane  M  (2  *  aMx  ♦  bMy  ♦  CM)  of  Figure  3-6  transforms  to  the 
point  (aM.hM,cM)  *n  abc  space,  and  planes  1,  J,  and  K.  transform  to  the 
points 

P|  *  (a|,b|,C|),  Pj  *  (aj.bj.Cj),  and  P^  *  (a^.b^.c^). 

The  inequality  (d)  (of  the  previous  section)  for  redundancy  of  half-space 
M  (point  (aM.bM.cM))  can  be  rewritten 

CM  *  (“px)aM  *  (“*VbM  +  p2* 

The  interpretation  in  abc  space  is  that  point  (»M,bM,cM)  lies  above  the 
plane 

c  ■  (-Px)a  *  (-Py)b  ♦  Pz. 

By  inspection  of  Equation  (3)  we  also  see  that  this  plane  is  determined  by 
the  three  points  Pj,  Pj,  and  P^.  The  three  “betweenness  of  slopes" 
conditions  (5)  define  three  vertical  planes,  each  of  which  passes  through 
two  of  tho  three  points  P|,  Pj,  and  PK.  Since  half-spaces  l,  J,  and  K  are 
all  UPPER  half-spaces,  the  set  of  redundant  points  PM  Is  a  bounded 


*0antsig  [28]  u-,e»  the  above  transform  in  the  context  of  linear  programming  and  Huffman  [53]  wo*  an  almost 
identical  transform  for  on  analysis  of  polyhedral  scenes.  Kanade  [57]  use*  this  transform  for  what  ho  calls  tho 
“origami  wotki*  of  tlveo -dimensional  figures. 
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region  of  nbc  space.  The  intersection  of  the  three  half-spaces  of 
conditions  (5)  and  the  half-space  (4)  is  therefore  the  region  in,  on,  or 
above  the  triangle  determined  by  points  P j,  Pj,  and  P^.  □ 

Corollary:  Let  the  N  UPPER  half-spaces  l^,  I Ig*  ...  11^  transform  to 
points  P-|,  P2,  .  .  .  P«j  by  Equation  (G).  A  half-space  H(  is  redundant  iff 
point  P|  does  not  lie  on  the  bottom  part  of  the  convex  hull  of  points  P^, 

p2-  •  •  •  PN- 

We  can  construct  the  (bottom  part  of  the)  convex  hull  of  N  three-dimensional 
points  in  0(N  log  N)  time.  We  do  this  by  augmenting  the  algorithm  of  Preparata  and 
Hong  [83]  for  constructing  the  (entire)  convex  hull  with  a  mechanism  for  separating 
the  bottom  part  of  the  hull  from  the  top.  One  way  to  do  this  is  to  maintain  with  each 
face  F  of  the  convex  hull  a  vector  perpendicular  to  F  that  points  toward  the  inside 
(as  opposed  to  the  outside)  of  the  hull.  The  bottom  faces  of  the  hull  are  those 
faces  whose  vectors  point  upward  and  the  vectors  for  the  top  faces  point 
downward.  (Note  that  there  will  be  some  vertices  in  both  the  top  and  bottom  parts 
of  the  hull.  These  are  the  vertices  that  bound  both  top  and  bottom  faces.)  The 
bottom  and  top  parts  of  the  convex  hull  are  therefore  separable  in  0(N)  time,  once 
the  entire  convex  hull  is  constructed  in  0(N  log  N)  time. 

— > 


Figure  3-7:  Transform  of  a  convex  hull. 


To  find  the  intersection  of  the  UPPER  half-spaces,  we  must  transform  the  bottom 
part  of  the  convex  hull  to  xyz  space.  The  0(N)  vertices  transform,  of  course,  to 
planes.  But  there  is  much  more  information  than  that  in  the  convex  hull.  For 
instance,  defining  each  face  of  the  hull  there  are  three  coplanar  vertices  I,  J,  and  K. 
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3  The  plane  that  these  vertices  define  is  transformed  to  a  point  in  xyz  space.  This 
point  is  where  pi, mcs  I,  J,  and  K  (of  xyz  space)  intersect.  (This  follows  from  the 
fact  that  the  transform  preserves  incidence  between  points  and  planes.)  Also,  as 
illustrated  in  figure  3-7,  if  faces  and  of  the  (bottom  part  of  the)  convex  hull 
share  an  edge  VjVj,  then  in  xyz  space  faces  V(  and  Vj  share  an  edge  F^FL.  In  fact, 
even  the  unbounded  faces  of  U  in  xyz  space  can  be  obtained  from  the  transform. 
These  faces  of  U  correspond  to  vertices  at  the  boundary  between  the  top  and 
bottom  parts  of  the  convex  hull  in  abc  space.  Little  computation  is  thus  required  to 
construct  U  after  the  transform  from  abc  space  since  all  the  faces,  vertices,  edges, 
etc.  are  directly  obtainable  from  the  transform.  We  have 

Theorem  7^  The  intersection  of  N  UPPER  three-dimensional  half-spaces 
can  be  constructed  in  0(N  log  N)  time. 

Proof:  Since  the  transform  costs  only  0(N)  time,  the  total  time  to 
construct  U  is  dominated  by  the  lime  to  construct  the  convex  hull  in  abc 
space,  which  is  0(N  log  N)  time.  □ 

3.1 .5.  Intersecting  Half-spaces  in  Four  or  More  Dimensions 

Suppose  that  in  K  dimensions  we  can  construct  the  convex  hull  of  N  points  in 
H(N,K)  time;  we  can  then  construct  the  intersection  of  N  UPPER  K-dimensional 
half-spaces  in  0(l!(N,K))  time.  The  algorithm  is  a  straightforward  extension  of  the 
one  that  we  used  in  the  three-dimensional  case.  We  first  transform  the  N 
half-spaces  to  N  points  in  K-space.  Then  we  construct  the  convex  hull  of  these  N 
points  in  H(N,K)  time.  Then  we  partition  the  top  from  the  bottom  part  of  the  hull  and 
transform  the  bottom  part  back  to  obtain  the  intersection  of  the  N  UPPER 
half-spaces,  of  the  hull,  and  transform  the  bottom  (top)  part  of  the  hull  back. 

In  this  section  we  present  an  algebraic  description  of  the  components  of  the 
K-dimensional  algorithm.  The  first  step  is  to  describe  algebraically  the  conditions 
under  which  an  UPPER  half-space  M  is  redundant  with  respect  to  a  set  S  of  K  UPPER 
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half-spaces.  We  then  present  a  general  duality  transform  for  K-space  and  interpret 
the  conditions  for  redundancy  in  the  transform  space.  Finally,  we  characterize 
redundancy  among  N  UPPER  half-spaces  in  terms  of  the  convex  hull  of  the  N  points 
to  which  they  transform. 

3.1. 5.1.  Algebraic  Description  of  Redundancy 

We  must  first  introduce  some  terminology.  Let  a  "j-face"  of  a  K-dimensional 
polytopc  be  denoted  nr.  follows:  a  vertex  is  a  0-face,  an  edge  (a  line  segment)  is  a 
1-face,  etc.  Let  S  be  a  set  of  K  UPPER  K-dimensional  half-spaces  defined  by 

K-1 

XK  *  £  Vl  +  i=1 . K 

j=1 

and  let  M  be  an  UPPER  K-dimensional  half-space  defined  by 

K-1 

XK  *  X  aM,*j  +  aMK- 
j*1 

Let  the  matrices  B  and  C  be  defined  by 


an  a,2  . 

•  au-i  *1 

"«1k" 

B  * 

a21  a22  * 

•  a2,K-1  *1 

and  C  =  - 

a2k 

_aK1  aK2  - aK,K-1 

_aKK_ 

The  flats  bounding  the  K  half-spaces  of  S  meet  at  a  point  Pa(P-|.P2 . Pfc) 

defined  by 

BP' =  C.  (7) 

In  three  dimensions  the  planes  I,  J,  and  K  determine  three  rays  rjj,  rj^,  and  r^j 
originating  at  point  P.  In  K  dimensions  we  have  K  such  rays,  denoted  r -j , 
where  rj  is  the  ray  determined  hy  the  K-1  elements  of  S  -  {i}.  If  (aj-j,  a\2<  •  •  •  Oj^) 
is  a  vector  pointing  in  the  same  direction  as  ray  i,  then  we  may  write  ray  i  in 
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parametric  form  as 

(x ] .Xg.  XK)  *  (Pi  *^2'  *  u  <XJK^*  **  0. 

On  the  other  hand,  each  point  (xi^.  ...  xK)  of  ray  i  must  also  lie  on  all  K  >  1  flats 
of  S  -  (i): 


Combining  these  two  equations  we  obtain 


To  solve  this  system  of  K  -  1  equations  in  K  unknowns  we  take  advantage  of  the 
following  properly  of  cofactors; 

K 

2  b;j  cofmj(B)  »  det(B),  if  i  «  m  (8) 

i*1  *  0  if  i  m. 

where  cofjj(B)  is  the  cofactor  of  b,j.  It  follows  that  the  general  solution  is 
(oc{1*ac{2'  *  *  **  ®iK^  *  ^  (cofji(B),  cof^CB). ....  cofjj^(B)), 
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where  is  a  constant  for  each  i.  The  choice  of  Cj  is  not  entirely  arbitrary  since 
SGN(Cj)  determines  whether  ray  i  points  up  or  down.  The  correct  choice  for  Cj 
satisfies 

K-1 

p«,+  a,K  <  2  a,j  (pj  4  a,j)  +  aiK 

j=1 

because  each  point  of  ray  i  (beyond  point  P)  should  lie  below  the  half-space  i. 
Using  the  definition  of  P  and  the  above  properties  of  cofactors  we  can  prove  that 
0  £  Cj  del(B).  The  vectors  a  are  therefore  chosen  to  be 

(«i1’«i2’  •  •  -«iK^  =  x  (cofji (B),  cofj2(B) . cofjK(B)). 

Theorem  Let  the  UPPER  half-space  M,  the  set  S  of  K  UPPER 
K-dimenoionai  half-spaces,  point  P,  and  K  rays  rj  be  as  given  above. 
Half-space  M  is  redundant  with  respect  to  the  K  half-spaces  of  S  iff 

K-1 

PK  *  I  aMjpj  4  aMK  <9> 

"  j=1 


and 


K-1 

X  <Vmj  -  «.*•  i=1 . (10) 

j=1 


Proof:  Half-space  M  is  redundant  iff  flat  M  lies  above  the  rays  rj,  i  = 
1 ,  .  .  K.  We  can  express  this  condition  as 


K-1 

PK  4  l'«,K  *  S  aMj  <Pj  4  »«„>  4  aMK’  V»  >  is1 . 

j=1 

If  M  is  redundant,  then  we  can  obtain  condition  (9)  by  simply  setting  u  = 
0.  Condition  (10)  results  from  dividing  by  u  and  taking  the  limit  as  u  -*  oo. 
Conversely,  if  conditions  (9)  and  (10)  are  satisfied,  then  M  must  be 
redundant.  Simply  multiply  condition  (10)  by  u  >  0  and  add  the  result  to 
condition  (9).  □. 
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3. 1.5.2.  The  General  Transform 
• 

Before  we  can  characterize  redundancy  in  the  transform  space,  we  must  first 
elucidate  the  important  properties  of  the  K-dimensional  transform.  This  transform 
maps  not  only  flats  to  points  and  points  to  flats,  but  also  j-spaces  to  K-j-1  -spaces. 

Gf!ner<l1  Transform:  The  transform  of  a  j-dimensional  subspace  of 
K-space  is  the  set  of  all  flats  that  contain  it.  This  is  a  K-j-1 -dimensional 
subspace  of  flats.  Since  each  flat  can  be  readily  represented  as  a  point, 
however,  the  transform  of  the  j-dimensional  subspace  is  represented  as  a 
K-j-1 -dimensional  subspace  of  points. 

Theorem  9:  The  general  transform  preserves  incidence.  In  other 
words,  if  a  j-j -dimensional  subspace  of  K-space  is  a  subspace  of  a 
j2-dimensional  subspace,  then  the  transform  of  the  j2~dimensional 
subspace  is  a  subspace  of  the  transform  of  the  j-j  -dimensional  subspace. 


Proof:  We  may  interpret  the  ji -dimensional  subspace  as  an 
intersection  of  K-j  -j  flats  and  the  j2*dimensional  subspace  as  an 
intersection  of  K-j2  flats.  That  is,  we  may  define  define  the 
j-j  -dimensional  subspace  by 


a11  al2 

•  •  • 

_  a1K 

a2l  ’  a22 

•  •  • 

a2.K-1  ’1 

a2K 

•  • 

• 

•  • 

X1 

• 

•  • 

• 

•  • 

X2 

• 

•  • 

• 

•  • 

. 

3  • 

• 

aK-i2.l  aK-j2.2 

•  •  • 

aK-J2,K-1  *1 

* 

aK-j2,K 

•  • 

•  * 

• 

• 

_XK_ 

• 

• 

_aK-j,.1  aK-j,,2 

*  *  * 

aK-it.K-1  *1_ 

_aK-jl.K_ 

and  the  j2*dimensional  subspace  by 
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The  transform  of  the  j-j -dimensional  subspace  is  the  set  of  points  of  the 
form 


and  the  transform  of  the  J2 -dimensional  subspace  is  the  set  of  points  of 
the  form 


It  is  easy  to  see  that  the  transform  of  the  j2-dimensional  subspace  is  a 
subspace  of  the  transform  of  the  j-j -dimensional  subspace.  □. 
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3. 1.5. 3.  Redundancy  in  the  Transform  Space 

Having  presented  both  a  characterization  of  redundancy  of  a  half-space  and  the 
general  duality  transform,  we  can  now  interpret  the  conditions  for  redundancy  in  the 
transform  space.  The  K  flats  bounding  the  half-spaces  of  S  transform  to  K  points  by 
the  formulas 


K-1 


j»1 


aUxi  +  a* 


(aa,  al2, 


3:1/),  i  ®  1  |  ■  ■  .,  K, 


and  the  flat  bounding  half-space  M  transforms  to  point  M  by 
K-i 

XK  *  S  aMjxj  +  aMk  **  ^aM1*  aM2’  *  *  •»  aMK^ 

J*1 

Since  the  transform  preserves  incidence  (Theorem  9),  the  flat  determined  by  the  K 

points  (fljj.ajg . aiK),  i  =  1, . . K  is  simply  the  transform  of  the  point  P 

determined  by  the  K  flats  of  S  (Equation  (7)).  Letting  the  coordinates  of  the 
transform  space  be  z-j ,  Z2, . . .  ,  z^,  this  flat  is 

K-1 

2k  s  Z  (-Pj)  ♦  Pk* 

•  jSl 

We  can  now  interpret  the  first  condition  (9)  for  redundancy  of  half-space  M  (with 
respect  to  S)  as  a  condition  on  point  M  and  the  flat  determined  by  the  transforms  of 
the  K  elements  of  S:* 


K-1 

aMK  *  X  (*P;)  aMj  +  PK' 

j*1 

We  have  therefore  proved 

Theorem  10:  The  first  condition  (9)  for  redundancy  of  UPPER 
half-space  M  with  respect  to  the  K  UPPER  half-spaces  of  S  requires  that 
point  M  (in  the  transform  space)  lie  above  the  flat  determined  by  the  K 
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points  (aj1t  ai2 . aiK),  i  =  1,  . .  K. 

The  other  conditions  for  redundancy  of  M  (the  "betweenness  of  slopes" 
conditions  (10))  map  to  a  set  of  K  half-spaces  in  K-1 -space 

K-1 

X  *ij2j  -  «iK 

i=1 

in  which  the  point  aM2<  •  •  •  <  aM,K-1^  must  lie*  Tlie  (K-1 -dimensional)  points 

that  satisfy  these  K  conditions  are 


"  2, " 

22 

a11  a21  • 
a12  a22  1 

•  aK-1,1 

•  aK-1,2 

"  U1  " 

U2 

• 

_alK  a2X  • 

•  aK-1,K_ 

where 

K 

X  u,  *  1  and  u,>  0,  i=1 . K.  (12) 

i=1 

That  is,  the  points  that  satisfy  the  "betweenness  of  slopes"  conditions  must  lie 
inside  the  convex  hull  of  the  K-1 -dimensional  points 

(aji,  ajg . ai,K-1^' '  =  1.  •  •  •«  K.  We  prove  this  assertion  as  follows.  Any  point  in 

K-1 -space  can  be  written  in  the  form  of  equation  (11)  if  the  restrictions  of  (12)  are 
ignored.  (We  have,  in  fact,  one  extra  degree  of  freedom.)  The  "betweenness  of 
slopes"  conditions  are  therefore 

K-1  K-1  K 

2  « ijZj  »  /L  det(B)  cofj (B)  *  <*;<.  s  det(B)  coftlc(B),  i=  1,...,K. 

j=1  j=1  h=1 


which  reduce  to 
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K  K-1 

dct(B)  X  uh  X  c°f, i.  dct(B)  cof!K(B),  i*1,...,|C 
h*1 *  j=1 

But  since  ail  entries  in  the  Kth  column  of  B  are  -1,  the  properties  of  cofactors  (8) 
allow  us  to  reduce  this  further  to 


*  K 

det(B)  cofjk(B)  2  uh  ♦  deUBJu, 
h®1 


>  det(B)  CfffIK(8),  i  =  1 . K 


(13) 


This  is  certainly  satisfied  for  all  i  If  the  restrictions  of  (12)  hold;  the  converse  is 
also  true.  Any  point  Q  that  does  not  satisfy  (12)  can  not  lie  in  the  convex  hull  of 

the  points  (a^,  aj2,  ....  ai,K-1  >* '  *  1 . K.  Here  we  can  take  advantage  of  the 

extra  degree  of  freedom  mentioned  above  to  express  the  point  Q  as  a  linear 
combination  of  these  K  points  such  that 

K  • 

I  u„  *  I  but  u,  <  0  for  some  i. 
h=1 

It  is  easy  to  see  that  such  a  point  does  not  satisfy  the  conditions  of  (13)  and  is 
therefore  not  redundant  with  respect  to  S. 


We  have  just  proved 

Theorem  1 1 ;  The  "betweenness  of  slopes"  conditions  (10)  for 
redundancy  of  UPPER  half-space  M  with  respect  to  the  K  UPPER 

half-spaces  of  S  require  that  the  point  (afvsi .  a^2 . aM,K- 1 )  ,ie 

the  convex  hull  of  the  points  (ajj,  ai2,  •  • ..  a]^),  •  *  1 .....  K.. 

Combining  Theorems  1 0  and  1 1  we  obtain 

Theorem  1 2:  UPPER  half-space  M  is  redundant  with  respect  to  the  K 
UPPER  half-spaces  of  S  iff  point  M  lies  directly  above  some  point  in  the 
convex  hull  of  the  points  (0^,  ai2 . a,^),  i  *  1, ....  K. 
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3. 1.5.4.  Redundancy  Among  PJ  Half-spaces 


In  this  section  we  characterize  redundancy  of  half-space  M  with  respect  to  a  set 
of  N  >  K  K-dimensionnl  half-spaces.  (If  N  is  less  than  K  then  there  can  not  be  any 
redundant  half-spaces  unless  some  are  parallel.) 

Theorem  1 3:  Let  T  he  a  set  of  N  UPPER  K-tlimensional  half-spaces. 
UPPER  haif-rspace  M  is  redundant  with  respect  to  the  half-spaces  of  T  iff 
M  is  redundant  with  respect  to  the  half-spaces  of  a  subset  S  of  T  that 
contains  exactly  K,  half-spaces. 

Proof:  It  is  clear  that  if  M  is  redundant  with  respect  to  a  subset  S  of  T, 
then  it  is  redundant  with  respect  to  T.  We  shall  now  establish  that  the 
converse  is  also  true.  Let  U  be  the  intersection  of  the  N  UPPER 
half-spaces  of  T.  U  is  a  convex  polytope  because  it  is  an  intersection  of 
half-spaces.  If  a  half-space  M  is  redundant  then  its  boundary  (flat  M) 
lies  completely  above  U.  Let  V  be  the  point  of  U  that  is  closest  to  flat 
M.  (If  the  closest  point  is  not  unique,  then  let  V  be  any  vertex  of  U  that  is 
in  the  set  of  closest  points.)  Since  U  is  a  convex  polytope,  V  is  (or  can 
be  chosen  to  be)  a  vertex  of  U.  Let  S  be  the  set  of  half-spaces  whose 
boundaries  meet  at  point  V.  If  the  half-spaces  of  T  are  in  general 
position,  then  S  contains  exactly  K  half-spaces.  Since  U  is  convex,  we 
can  travel  from  vertex  V  along  the  boundary  of  U  in  any  direction  and  the 
distance  to  flat  M  will  be  nondecrcasing.  It  follows  that  the  boundary  of 
U  can  never  intersect  flat  M  and  therefore  half-space  M  is  redundant 
with  respect  to  the  K  half-spaces  of  S.  □. 

We  have  just  characterized  redundancy  of  a  half-space  M  with  respect  to  N 
half-spaces  in  terms  of  redundancy  with  respect  to  K.  half-spaces.  But  M  is 
redundant  with  respect  to  a  set  of  K  half-spaces  iff  point  M  lies  above  a  point  in 
the  convex  hull  of  the  points  to  which  the  K  half-spaces  transform.  We  have 
therefore 

Theorem  Id:  An  UPPER  half-space  M  is  redundant  with  respect  to  a  set 
T  of  N  UPPER  K-dimensional  half-spaces  iff  the  point  M  (to  which 
half-space  M  transforms)  lies  directly  above  a  point  in  the  convex  hull  of 
the  N  points  of  the  transform  of  the  N  half-spaces  of  T. 

Since  by  Theorem  9  the  components  of  the  convex  hull  correspond  to  components  of 

the  intersection  of  half-spaces,  we  have 
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Theorem  1 5:  If  H(N.K)  is  the  time  lo  construct  the  convex  hull  of  N 
points  in  K-space,  then  N  UPPER  K-dimensional  half-spaces  can  be 
intersected  in  0(H(N,K))  time. 


3.1.6.  Open  Problems 


There  ore  still  a  few  open  problems  concerning  intersection  of  half-spaces. 

1.  The  time  to  intersect  N  (UPPER)  K-dimensional  half-spaces  depends  on 
the  time  ll(N,K)  to  construct  the  convex  hull  of  N  K-dimensional  points. 
We  know  that  the  worst-case  time  complexity  is  0(N  log  N)  for  two 
and  three  dimensions  but  for  four  or  more  dimensions  only  the  Q(N^) 
time  lower  bound  has  been  proven.  (See  Section  1.1.1  for  a 
description  of  our  knowledge  of  convex  hulls.)  For  K  >  4  dimensions 
we  lack  tight  upper  and  lower  bounds  on  H(N,K).  ' 

2.  Under  some  conditions  the  expected-time  for  intersection  of 
half-spaces  may  be  less  than  the  worst-case  time.  We  may  use  fast 
expected-time  convex  hull  algorithms  to  obtain  fast  expected-time 
intersection  of  half-space  algorithms.  For  example,  in  two  and  three 
dimensions  if  the  expected  number  of  nonredundant  half-spaces  is 
0(Mp)  for  some  p  <  1,  then  N  half-spaces  can  be  intersected  in  O(N) 
expected-time  [16].  Since  the  convex  hull  of  N  K-dimensional  points 
can  be  constructed  in  0(N)  expected-time,  if  the  K  coordinates  have 
independent  distributions  [10,  30]  ,  we  can  intersect  N  K-dimensional 
half-spaces  in  O(N)  expected-time  if  the  N  half-spaces  transform  to  N 
points  whose  K  coordinates  are  distributed  independently.  Under  what 
other  conditions  may  we  intersect  half-spaces  in  fast  expected-time? 

3.  We  can  intersect  two  convex  (three-dimensional)  polyhedra  in  0(N  log 
N)  time  by  simply  treating  it  as  a  problem  of  intersecting  half-spaces. 
Can  we  improve  this  to  0(N)  time?  (In  two  dimensions,  convex  N-gons 
can  be  intersected  in  0(N)  time  whereas  it  requires  0(N  log  N)  time  to 
intersect  N  half-planes  in  the  worst  case.) 

4.  How  fast  can  we  intersect  N  half-spaces  on-line  in  two  or  more 
dimensions?  (Shamos  [91]  has  presented  an  0(N  log  N)  time  planar 
on-line  convex  hull  algorithm  and  Preparata  [78]  has  refined  that  to  an 
0(N  log  N)  time  real-time  algorithm.  Qoth  of  these  algorithms  update 
the  convex  hull  as  each  point  is  read  —  rather  than  operating  on  all  N 
points  collectively  —  but  the  on-line  algorithm  may  require  up  to  0(N) 
time  for  any  particular  update  whereas  the  real-time  algorithm  always 
requires  at  most  0(log  N)  update  time.) 
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3.2.  Union  and  Intersection  of  Disks 

We  present  the  problem  of  constructing  the  union  or  intersection  of  a  set  of  disks 
(the  interiors  of  a  set  of  circles)  not  for  its  applications,  but  because  the  solution 
presents  several  techniques  that  are  useful  for  solving  many  other  problems.  We 
introduce  two  new  transforms,  inversion  and  embedding  in  a  higher  dimension. 
Inversion  is  used  to  convert  problems  that  involve  circles  or  spheres  to  problems 
that  involve  lines  or  planes.  Embedding  in  a  higher  dimension  adds  another  degree 
of  freedom  to  the  problem,  which  can  permit  application  of  techniques  not  applicable 
to  the  original  problem.  The  disk  algorithm  combines  inversion  and  embedding  to 
transform  the  problem  of  constructing  the  union  or  intersection  of  a  set  of  N  disks  to 
the  problem  of  intersecting  N  three-dimensional  half-spaces,  which  we  have  solved 
in  0(N  log  N)  time.1®  Another  important  feature  of  the  algorithm  is  that  it 
demonstrates  an  example  of  how  a  nonconvex  object  (the  union  of  disks)  can  be 
represented  by  a  convex  object  (an  intersection  of  half-spaces). 

3.2.1.  Representation  of  Disks  and  Their  Union  or  Intersection 

A  disk  is  the  set  of  all  points  within  a  given  positive  radius  R  from  a  planar  point  P. 
If  the  coordinates  of  P  are  (Px,Py),  then  the  disk  may  be  represented  as  a  triple 
(Px,Py,R)  and  a  set  of  disks  as  a  set  of  such  triples.  The  best  representation  for 
the  union  or  the  intersection  of  a  set  of  disks  depends  on  what  kinds  of  information 
about  the  union  or  intersection  we  want  to  retrieve  efficiently.  For  several 
applications  the  best  representation  will  be  as  an  intersection  of  half-spaces  in  the 
transform  space  (to  be  described),  but  for  other  applications  we  may  need  a 
representation  in  the  original  space.  We  will  not  describe  either  representation  in 
detail  because  we  have  already  described  the  representation  for  intersection  of 
half-spaces  and  the  representation  in  the  original  space  is  only  a  slight  modification 


®Shamos  and  Hooy  [9/|]  describe  an  0(N  log  N)  time  intersection  ol  N  disks  that  is  a  modification  to  their 
algorithm  for  intersecting  half -planes.  Their  algorithm  unfortunately  does  not  extend  to  an  0(N  log  N)  time 
algorithm  for  the  union  of  N  disks. 
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of  the  representation  for  polygons. 


3.2.2.  Lower  Bound  for  the  Union  or  Intersection  of  Disks 


We  prove  that  an  algorithm  that  constructs  the  union  or  intersection  of  N  disks 

can  be  used  to  sort.  Since  there  is  an  Q(N  log  N)  time  lower  bound  for  sorting 

(under  the  same  model  of  computation  used  in  the  construction  below),  we  have  an 

Q(N  log  N)  time  lower  bound  for  constructing  the  union  or  intersection  of  N  disks. 

Theorem  16:  The  construction  of  the  union  or  the  intersection  of  N 
disks  takes  ft(N  leg  N)  time. 


h 

o 


Proof:  Let  S  be  a  set  of  N  reals  tj  *  [0,1),  i  a  1 . N.  For  each  real  tj 

we  have  a  disk  dj  centered  at  (xj.y;)  with  radius  two  where 

(Xj,  yj)  »  (  cos(2ntj),  sin(2TTtj)). 


As  illustrated  in  Figure  3-8,  the  points  (Xj.yp  are  all  on  the  unit  circle  and 
the  boundaries  of  the  disks  d-,  are  tangent  to  the  unit  circle.  The  union  of 
the  disks  dj  is  represented  by  a  closed  chain  of  circular  arcs 
(ab-bc-cd-da.in  Figure  3-8).  The  order  of  the  arcs  in  this  chain  forms  a 
sort  for  tiie  N  reals  t;.  Similarly,  an  intersection  of  the  same  N  disks  is 
represented  by  a  closed  chain  of  circular  arcs  (he-ef-fg-gh  in  Figure  3-8) 
whose  order  sorts  the  reals  tf. 


Figure  3-8:  Sorting  with  a  union  or  intersection  of  disks  algorithm. 
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The  proof  of  the  lower  bound  requires  an  Q(N  log  N)  time  lower  bound 
for  sorting  under  the  same  model  of  computation  used  to  construct  the  N 
disks  dj  from  the  M  reals  tj.  Since  the  construction  of  circles  uses  the 
functions  sine  and  cosine,  a  model  of  computation  with  only  linear 
functions  of  the  input  is  not  sufficient.  As  for  the  intersection  of 
half-spaces,  we  may  use  Friedman's  [43]  result  that  sorting  has  an 
Q(N  log  N)  time  lower  bound  even  when  arbitrary  functions  are  allowed  at 
internal  nodes  of  the  decision  tree  and  the  output  functions  are  analytic. 
Since  Sine  and  cosine  are  analytic,  the  Q(N  log  N)  lower  bound  for  sorting 
carries  over  to  the  union  and  intersection  of  disks. 


3.2.3.  The  Inversion  Transform 

Our  algorithm  for  constructing  the  intersection  or  union  of  N  disks  is  based  on  the 
properties  of  the  inversion  transform.  We  will  first  describe  the  two-dimensional 
transform  and  then  generalise  to  three  (and  higher)  dimensions.  For  more  information 
on  inversion  we  refer  the  reader  to  Dodge  [36]. 

The  inversion  transform  is  determined  by  two  parameters:  (1)  the  center  of 
inversion,  and  (2)  the  radius  of  inversion.  For  simplicity  of  exposition  we  shall 
assume  (for  now)  that  the  center  of  inversion  is  the  origin  and  that  the  radius  of 
inversion  is  one.  If  a  point  P  has  polar  coordinates  (R,0),  then  the  inversion 
transform  of  P  is 

(  R,  6  )  -*■  (  1/R.0  ). 

Inversion  maps  a  vecto:  in  the  direction  Q  to  another  vector  in  the  same  direction 
but  with  its  magnitude  "inverted."  Note  that  inversion  is  involutory  —  application  of 
inversion  twice  yields  the  original  point.  Figure  3-9  illustrates  another  important 
property  of  inversion  in  the  plane.  A  circle  that  passes  through  the  center  of 
inversion  transforms  to  a  line  that  docs  not  pass  through  the  center  of  inversion,  and 
vice  versa.  Furthermore,  the  interior  of  the  circle  transforms  to  one  of  the 
half-planes  determined  by  that  line  and  the  exterior  of  the  circle  transforms  to  the 
other  half-plane.  The  properties  of  inversion  in  three  dimensions  are  analogous.  The 
transform  can  be  expressed  in  spherical  coordinates  as 

(  R,  0,  $  )  ■+  (  1/R,0,$  ). 
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This  transform  is  involutory,  as  in  two  dimensions,  and  it  also  transforms  any  sphere 
that  passes  through  the  center  of  inversion  to  a  piano  that  does  not  pass  through 
the  center  of  inversion.  The  interior  of  the  sphere  transforms  to  a  half-space 
bounded  by  that  plane  and  the  exterior  of  the  sphere  transforms  to  the  other 
half-space. 


Figure  3-9:  Inversion  transforms  lines  to  circles  and  vice  versa. 


3.2.4.  Algorithm  for  Intersection  or  Union  of  Disks 

We  will  now  exploit  the  properties  of  the  inversion  transform  to  construct  a  fast 
algorithm  for  the-  union  or  intersection  of  disks.  We  will  first  describe  the  simple 
case  where  the  circles  bounding  the  N  disks  share  a  common  point  P,  and  then 
generalize  the  result  to  an  arbitrary  set  of  N  disks  in  the  plane. 

Figure  3-10  illustrates  the  special  case  where  the  N  boundary  circles  share  a 
common  point  P.  Let  Cj  denote  disk  i,  for  i»1  ,...N.  Since  circle  i  passes  through  point 
P.  inversion  about  P  transforms  Cj  to  o  half-plane  Hj.  It  follows  that  the  union  of  the 
N  disks  transforms  to  the  union  of  N  hnif-planes: 

N  N 

U  Cs  -  U  H;. 
i«1  i»1 


Similarly,  if  X  denotes  the  complement  of  X,  then  the  union  of  the  Cj  transforms  by 
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N  N  N  ~ 

u  c  -  u  h,  *  n  c;. 

i=1  i=1  i=1 

which  can  also  be  solved  as  an  intersection  of  half-planes.  Since  N  half-planes  can 

be  intersected  in  0(N  log  N)  time,  we  have 

Theorem  ±7:  The  union  or  intersection  of  N  disks  can  be  constructed  in 
0(N  log  N)  time  when  the  circles  bounding  the  N  disks  share  a  common 
point  P. 

In  general,  however,  the  N  circles  will  not  have  any  point  P  in  common.  In  fact, 
the  disks  may  he  entirely  disjoint.  Nevertheless,  we  can  manipulate  the  general 
problem  so  that  it  looks  sufficiently  like  the  special  case  that  inversion  can  be 
applied  to  obtain  a  fast  algorithm. 


Figure  3-10:  Special  case  for  inversion:  AU  boundary  circles  meet  at  point  P. 

Theorem  1 8:  The  union  or  intersection  of  N  disks  can  be  represented 
as  a  convex  polyhedron  in  0(N  log  N)  time. 

Proof:  We  illustrate  the  construction  in  Figure  3-11.  We  first  embed 
the  N  disks  in  three  dimensions  with  the  disks  all  located  in  the  xy  plane. 
We  then  choose  an  arbitrary  point  P  that  does  not  lie  in  the  xy  plane.  For 
each  disk  c  there  is  a  unique  sphere  that  posses  through  point  P  and  that 
intersects  the  xy  plane  at  circle  c.1  ^  We  can  thus  represent  the  N  disks 


^  A  sphere  is  determined  by  four  parameters,  the  three  coordinates  for  its  center  and  the  radius.  A  circle  and  a 
noncoptaror  point  determine  a  unique  sphere  because  requiring  the  sphere  to  pass  through  the  circle  costs  three 
degrees  of  freedom  and  requiring  the  sphere  to  pass  through  the  point  determines  the  fourth. 
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in  the  xy  plane  by  N  balls  whose  (spherical)  boundaries  share  a  common 
point  P.  Inversion  about, point  P  transforms  the  M  spheres  to  N  planes,  the 
balls  to  half-spaces,  and  and  the  exteriors  of  the  balls  to  complementary 
half-spaces.  The  intersection  of  M  disks  is  thus  represented  by  the 
intersection  of  N  half-spaces.  Similarly,  we  represent  the  union  of  N  disks 
by  (the  complement  of)  the  intersection  of  N  (complementary) 
half-spaces.  Since  we  can  intersect  N  half-spaces  in  0(N  log  N)  time,  we 
can  represent  the  union  or  intersection  of  N  arbitrary  planar  disks  by  a 
convex  polyhedron  in  0(N  log  N)  time.  □. 


Figure  3-11:  General  case  for  intersection  or  union  of  N  disks. 


3.2.5.  Related  problems 

We  will  now  briefly  describe  a  few  issues  related  to  the  union  of  disks  that  we 
have  not  yet  covered. 

1 .  We  con  generalize  our  technique  for  representing  the  union  or 
intersection  of  disks  to  also  allow  subtracting  circular  regions  from  a 
union  or  intersection  of  disks.  Since  each  circle  mop3  to  a  plane  in 
three-space,  there  are  two  three-dimensional  half-spaces  that  we 
may  associate  with  a  given  circle.  One  half-space  corresponds  to  the 
interior  of  the  circle  and  the  other  corresponds  to  the  exterior.  We 
may  thus  represent  any  intersection  of  the  interiors  of  circles  or  their 
complements  by  an  intersection  of  half-spaces. 

2.  Suppose  that  we  wanted  to  preprocess  N  disks  so  that  given  an 
arbitrary  planar  point  we  could  quickly  determine  if  the  point  lies  in  any 
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of  the  N  disks.  If  we  represent  the  union  of  the  disks  ns  an 
intersection  of  half-spaces,  then  this  problem  becomes  the  problem  of 
determining  if  a  point  in  three-space  lies  within  a  convex  polyhedron. 
This  problem  is  equivalent  to  locating  a  point  within  two  planar 
straight-line  graphs.^  As  wo  mentioned  in  the  introduction,  location 
of  a  point  within  a  planar  graph  of  size  N  can  be  done  in  0(log  N)  time 
given  0(N  log  N)  preprocessing  time  and  0(N)  storage  [70].  We  must 
also  note  that  location  of  a  point  in  a  set  of  circles  is  a  decomposable 
searching  problem  [5]. 


3.3.  Derivation  of  the  Point  /  Flat  Duality 

The  intersection  of  half-spaces  involves  a  point  /  flat  duality  and  convex  hulls 
whereas  the  union  of  disks  depends  on  inversion  (and  embedding  in  a  higher 
dimension).  In  this  section  we  show  that  these  techniques  are  closely  related  by 
deriving  the  point  /  flat  duality  as  a  limiting  case  of  inversion,  linear  transforms,  and 
a  circle  /  point  duality.  In  the  construction  we  demonstrate  the  relationship 
between  convex  hulls  and  the  union  of  disks  or  half-spaces. 


Figure  3-12:  The  union  of  (the  interiors  of)  circles  that  meet  at  a  point  P. 


In  the  previous  section  we  constructed  the  union  of  (the  interiors  of)  N  circles 


12 

Simply  break  the  convey  polyhedron  m;o  UPPER  ar-l  LOWER  p.vti  and  project  both  parts  orthographically  to 
planar  graphs  in  the  x.y  plane.  If  a  three -dirne'i-sional  poml  P  projects  lo  a  planar  point  in  region  R  of  the  UPPER 
graph  arid  region  S  of  the  LOWER  giaph,  then  P  lies  within  the  convex  polyhedron  iff  P  lies  below  face  R  and 
above  face  S  of  Ihe  polyhedron. 
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with  a  point  P  in  common  by  first  using  inversion  to  transform  the  interiors  of  the 
circles  to  half-planes.  In  Figure  3-12  we  demonstrate  an  alternate  solution: 

Construction  o'  the  Union  of  M  Disks  whose  Boundaries  Share  a  Point  P 

1.  Let  Cj.  i  *  1,  ....  N  be  a  set  of  disks  whose  boundary  circles  meet  at 
point  P.  Transform  each  of  the  N  disks  Cj  to  the  points  Pj  diametrically 
opposite  point  P.  (This  is  the  Circle  /  Point  duality.) 

2.  Construct  the  convex  hull  of  the  set  of  the  points  Pj,  i  »  1 . N 

augmented  with  point  P. 

3.  A  disk  Cj  is  nonredundant  in  the  union  of  the  N  disks  iff  its  transform  Pj 
is  a  vertex  of  the  convex  hull.  To  obtain  the  circular  arcs  defining  the 
union  we  simply  transform  the  vertices  of  the  convex  hull  back  to 
disks,  obtaining  the  endpoints  of  the  arcs  from  the  points  where 
neighboring  disks  intersect. 

Since  the  most  expensive  step  in  the  algorithm  is  the  construction  of  the  convex 
hull,  the  total  lime  is  0(N  log  N).  We  v/ill  omit  the  proof  that  the  algorithm  correctly 
constructs  the  union  of  the  N  circles  because,  as  we  shall  see,  it  is  essentially 
equivalent  to  our  proof  of  the  algorithm  for  intersection  of  half-spaces. 


- > 


Figure  3-13:  The  union  of  half-planes  transforms  to  the  union  of  disks. 


In  Figure  3-13  we  illustrate  the  union  of  a  set  of  half-planes,  each  of  which  does 
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not  contain  point  P.1^  Dy  invar  tiny  the  half-pianos  about  point  P  we  transform  the 
union  of  half-planer,  problem  to  a  union  of  disks  problem  where  the  circles  bounding 
the  disks  all  meet  at  point  P.  Cut  as  we  just  saw,  we  can  construct  the  union  of  the 
disks  by  transforming  them  to  points  and  constructing  the  convex  hull  of  those 
points  (and  point  P).  We  will  now  describe  the  transform  algebraically. 

To  make  the  mathematics  simpler  we  can  assume  that  the  point  P  lies  on  the 
negative  y-axis: 

p  =  (o,r),  n  <  o. 

Our  algorithm  inverts  each  half-plane  about  point  P,  obtaining  the  interior  of  a  circle, 
and  then  maps  the  circle  to  the  point  diametrically  opposite  P.  We  express  this  as 

y  =  ajX  +  bj  ->  (  -aj/(bj-R),  R  +  1  /(bj-R)  ),  i  =  1 . N. 

To  show  how  this  relates  to  the  point  /  fiat  duality  that  we  used  in  Section  3.1  for 
intersecting  half-planes,  wc  must  now  use  a  trick. 


Recall  from  Step  3  of  our  algorithm  for  constructing  the  union  of  the  disks,  a  disk 
Cj  is  nonrodundant  iff  its  transform  Pj  is  a  vertex  of  the  convex  hull.  If  the  convex 
hull  is  translated,  stretched,  or  rotated  the  vertices  will  still  remain  vertices  and  the 
points  inside  the  hull  will  remain  inside  the  hull.  In  particular,  the  transform 


fx'-j  TR  0  M  +  f  0  " 

Ly'l  "  [o  -r2J  L/J  +  Lr3-R. 


applied  to  the  point  P  and  points  Pj.  i  *  1 . N  will  not  affect  our  determination  of 

which  points  lie  on  the  convex  hull  and  which  do  not.  Our  new  transform,  inversion 
of  a  half-plane  about  point  P  followed  by  the  circle  /  point  duality  and  our  linear 
transform,  can  be  expressed  as 


15 

Wo  can  alternately  think,  of  the  union  of  half-planes  as  Ihe  complement  of  the  intersection  of  the 
complementary  half-planes. 
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y  =  a,x  +  b, 


-a,R  -b;R  ' 
b;  -  R  ’  b,  R  ‘ 


Note  also  that  point  P  transforms  by 

P  *  (O.R)  ■+  (0,  -R). 

Taking  the  limit  as  R  ■*  -w>  we  obtain 

y  =  ajx  +  bj  -♦  (aj.bj) 

and  for  point  P 

P  =  (O.R)  -*  (0,  too). 

The  vertices  of  the  convex  hull  of  tlie  points  (aj.bj),  i  s  1 . N  augmented  with 

(0,+oo)  are  the  vertices  on  the  bottom  port  of  the  convex  hull  of  the  points  (aj.bj) 

(and  (0,+(n)).  This  is  exactly  the  characterization  given  for  intersection  (or  union) 

of  half-planes  in  Section  3.1 .  We  have 

Theorem  1 9:  The  point  /  flat  duality  is  a  limiting  case  of  inversion 
followed  by  the  circle  /  point  duality  and  a  linear  transform. 


3.4.  Summary 

In  this  chapter  we  have  presented  several  important  transforms  and  techniques: 

-  A  point  /  flat  duality  — 

This  maps  points  to  flats  and  flats  to  points.  Since  there 
are  already  a  number  of  algorithms  for  point  problems,  this 
transform  finds  greatest  use  in  transforming  problems  that 
are  expressed  (or  expressible)  in  terms  of  flats  to  problems 
that  involve  points.  Two  of  the  important  properties  of  the 
duality  that  we  described  are 

(a)  it  preserves  above/beiowness  between  points  and  flats, 
and 

(b)  it  preserves  distance  between  points  and  flats  in  the  x^ 
coordinate,  and  thus  preserves  incidence. 
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-  Embedding  into  a  higher  dimension  ~ 

This  gives  another  degree  of  freedom  to  the  problem  that 
allows  application  of  techniques  not  applicable  to  the 
original  problem  --  circles  con  become  spheres,  lines  can 
become  planes,  etc.  Since  the  higher-dimensional  object 
has  a  degree  of  freedom  that  v/e  can  choose  arbitrarily,  the 
objects  can  be  chosen  to  conform  to  an  expecially  simple 
case  (all  N  spheres  have  a  point  in  common).  In  most  of  the 
applications  of  embedding  in  a  higher  dimension  (in  this 
thesis)  the  problem  is  first  expressed  in  terms  of  circles 
and  then  embedded  into  a  higher  dimension  and 
re-expressed  in  terms  of  spheres. 

-  Inversion  — 

Inversion  is  a  circular  transform;  circles  map  to  circles 
(where  a  line  is  considered  to  be  a  circle  of  infinite  radius). 

In  particular, 

(a)  a  circle  that  passes  through  the  center  of  inversion  maps 
to  a  line  that  does  not  pass  through  the  center  of 
inversion,  and 

(b)  the  interior  of  a  circle  that  contains  the  center  of  inversion 
maps  to  the  exterior  of  a  circle  that  contains  the  center  of 
inversion. 

In  three  dimensions  we  have  the  same  relationships 
between  spheres  and  planes  and  in  K  dimensions  between 
K-spheres  and  flats.  Inversion  is  also  involutory  — 
application  of  it  twice  yields  the  original  object.  The  main 
use  for  inversion  (in  this  thesis)  is  transforming  problems 
that  are  expressed  (or  expressible)  in  terms  of  circles  or 
spheres  to  problems  that  involve  lines  or  planes  and  for 
which  fast  algorithms  are  known. 

In  this  chapter  we  have  also  seen  an  application  of  convex  hulls  in  the  intersection 
of  half-spaces  and  in  the  relationship  between  inversion,  linear  transforms,  and  the 
point  /  flat  duality.  In  later  chapters  convex  hulls  will  be  applied  to  several  other 
problems  that  involve  a  "network"  of  linear  parts. 
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4.  Construction  of  Nearest  and  Farthest  Point  Diagrams 

Geometric  transforms  arc  important  tools  in  the  construction  of  many  different 
kinds  of  nearest  and  farthest  point  diagrams.  These  diagrams  include  (nearest  and 
farthest  point)  Vornnoi  diagrams  in  Euclidean  and  spherical  spaces,  and  the  (nearest 
and  farthest)  edge  diagrams  of  a  convex  polygon.  Each  of  these  diagrams  is  a 
tesselation  of  space  into  sets  of  points  closest  to  (or  farthest  from)  the  elements 
(points  or  edges)  defining  the  diagram.  We  will  find  in  each  case  that  it  is  useful 
first  to  express  the. problem  in  terms  of  a  set  of  circles  that  define  the  diagram  and 
then  to  embed  the  problem  into  a  higher-dimensional  Euclidean  space,  which  allows 
us  to  use  techniques  not  applicable  to  the  original  problem.  In  the  following  sections 
we  describe  tiiese  diagrams  and  algorithms  for  constructing  them. 

4.1.  Euclidean  Voronoi  and  Delaunay  Diagrams 

Voronoi  diagrams  (also  catted  Thicsscn  diagrams  or  Oirichlet  tesselations)  find 
application  in  cluster  analysis  [51],  construction  of  contour  maps  [29],  construction 
of  Euclidean  minimal  spanning  trees  [93],  crystal  growth  [46],  and  several 
interesting  problems  in  geometry  [91].  Wc  can  easily  show  an  Q(N  log  N)  time 
worst-case  lower  bound  by  demonstrating  that  any  algorithm  that  constructs  a 
Voronoi  diagram  can  be  used  to  sort  [91];  the  challenge  is  to  construct  an 
0(N  log  N)  lime  algorithm.  Shainos  [89]  describes  on  0(N  log  N)  time 

divide-and-conquor  algorithm  for  construction  of  the  planar  Euclidean  Voronoi 
dianram  and  I  ee  and  Wo  no  F671  describe  an  OfN  loo  N)  time  aloorithm  for  the  I  «  and 

V  Vi  J  •  V  '  V  | 

Lf0  metrics  in  the  plane.  Orysdnle  and  Lee  [37]  present  an  0(N  c^°9  time 
algorithm  for  the  Voronoi  diagram  of  line  segments  (and  other  planar  objects),  which 
they  have  improved  to  0(N  log?N)  time.  Kirkpatrick  [59]  presents  an  0(N  log  N)  time 
algorithm  for  constructing  the  Voronoi  diagram  of  N  planar  line  segments.  Shamos 
[89],  Lee  and  Preparata  [66],  Prcparata  [80],  and  Lipton  and  Tarjan  [70]  have 
produced  fast  algorithms  for  searching  a  Voronoi  diagram  (or  any  other  straight-line 
planar  graph). 

The  algorithm  for  construction  of  a  Euclidean  Voronoi  diagram  that  we  describe 
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below  [23]  is  not  only  a  very  useful  result  in  itself,  but  it  also  serves  as  an  example 
of  the  use  of  several  important  algorithmic  tools.  We  use  the  technique  of 
embedding  into  a  higher  dimension  and  applying  inversion  (as  in  the  algorithm  for 
union  of  circles)  and  we  also  use  a  convex  hull  algorithm  (as  in  the  intersection  of 
half-spaces). 

4.1.1.  Definition  of  Planar  Voronoi  and  Delaunay  Diagrams 

Let  S  be  a  set  of  N  planar  points  such  that  no  four  points  are  cocircular.^  A 
nearest  point  planar  Voronoi  diagram  of  S,  as  pictured  in  Figure  4-1,  is  a  polygonal 
network  of  N  regions.  For  each  point  Pj  of  S,  region  Rj  is  the  set  of  all  points  of  the 
plane  that  are  closer  to  point  Pj  than  to  any  of  the  other  N-1  points  of  S.  Given  an 
arbitrary  point  P  in  the  plane,  we  can  Huts  determine  which  of  the  N  points  of  S  is 
closest  to  P  by  determining  which  of  the  M  regions  contains  point  P.  The  vertices  of 
these  polygonal  regions  are  called  Voronoi  points  and  the  polygonal  boundaries  of 
the  regions  are  called  Voronoi  polygons.  If  a  Voronoi  polygon  is  bounded,  then  it  is 
constructed  entirely  from  edges  of  the  Voronoi  diagram.  If  it  is  unbounded,  then  it 
includes  two  rays  of  the  diagram. 

Each  Voronoi  point  V  of  the  nearest  point  diagram  is  equidistant  from  the  three 
points  of  S  that  are  nearest  V.  This  yields  a  property  of  Voronoi  diagrams  that  we 
exploit  in  the  algorithm  of  Section  4.1.3. 

A  point  V  is  a  Voronoi  point  (of  the  nearest  point  Voronoi  diagram  of  S) 
iff  it  is  the  center  of  a  circle  that  passes  through  three  points  of  S  but 
does  not  contain  any  of  the  other  N  -  3  points  of  S. 

The  edges  of  a  Voronoi  diagram  connect  pairs  of  Voronoi  points  whose  corresponding 
circles  meet  at  two  common  points  of  S.  The  rays  are  determined  similarly  by  one 
Voronoi  point  from  the  nearest  point  diagram  and  one  from  the  farthest-point  Voronoi 

diagram  (described  below). 

I 

! 


^Sineo  nny  three  noncollincar  planar  points  determine  a  unique  circle,  four  planar  points  are  eocircular  only  in 
degenerate  cases. 


Figure  4-2:  F 


it  Point  Voronoi  Diagram 
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regions,  but  region  Rj  is  the  sot  of  all  points  in  the  plane  that  are  further  from  point 
Pj  than  any  other  point  of  S.  As  for  the  nearest  point  diagram,  there  is  a  set  of 
circles  that  define  the  Voronoi  points  for  a  farthest  point  diagram. 

A  point  V  is  a  Voronoi  point  of  the  farthest  point  Voronoi  diagram  iff  it  is 
the  center  of  a  circle  that  passes  through  three  of  the  points  of  S  and 
contains  nil  of  the  other  N  -  3  points. 

It  is  important  to  note  that  only  points  that  are  vertices  of  the  convex  hull  of  the  N 
points  of  S  have  nonempty  farthest  point  regions.  This  is  because  each  Voronoi 
point  V  of  the  farthest  point  diagram  must  be  equidistant  from  the  three  points  of  S 
that  are  farthest  from  V.  It  is  not  possible  to  construct  such  a  Voronoi  point  from 
points  of  S  that  are  not  on  the  convex  hull. 

\  I 

\  I 

\  I 

\  / 


Figure  4-3:  Planar  Delaunay  Diagram 

Both  the  nearest  and  farthest  point  Voronoi  diagrams  have  planar  straight-line 
duals,  called  Delaunay  diagrams.  Figure  4-3  illustrates  the  dual  of  a  nearest  point 
Voronoi  diagram  and  Figure  4-4  illustrates  the  dual  of  a  farthest  point  Voronoi 
diagram.  Both  of  these  dual  diagrams  form  a  triangulation  of  the  points  of  S.  The 
vertices  of  each  triangle  of  the  (nearest  point)  Delaunay  diagram  (Figure  4-3) 
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determine  a  circle  Hint  does  not  contain  any  of  the  other  N  -  3  points  of  S.  Similarly, 
the  vertices  of  each  triangle  of  the  dual  of  the  farthest  point  Voronoi  diagram 
(Figure  4-4)  determine  a  circle  that  contains  all  of  the  other  N  -  3  points  of  S. 


Figure  4-4:  Dual  of  a  Farthest  Point  Voronoi  Diagram 

Since  the  nearest  and  farthest  point  Voronoi  diagrams  and  their  duals  are  planar 
graphs,  the  number  of  Voronoi  points  is  at  most  2N  -  4  and  the  number  of  edges  is  at 
most  3N  -  G  for  N  >  2  [50].  Shamos  [39,  91]  and  Shamos  and  lloey  [93]  give  more 
information  on  Voronoi  and  Delaunay  diagrams. 

4.1 .2.  Representation  of  Voronoi  and  Delaunay  Diagrams 

The  representation  of  Voronoi  and  Delaunay  diagrams  should  enable  us  to  access 
conveniently  all  of  the  proximity  information  stored  in  the  diagrams.  This  does  not 
require  an  exotic  data  structure  —  we  can  accomplish  it  with  a  well-chosen  set  of 
arrays.  Since  the  rays  of  both  the  nearest  and  farthest  point  Voronoi  diagrams  are 
determined  by  one  nearest  Voronoi  point  and  one  farthest  Voronoi  point,  It  is  easiest 
for  us  to  describe  a  representation  for  both  the  nearest  and  farthest  point  diagrams 
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simultaneously.  One  such  representation  is  the  five-tuple 

(S,  V,  E.  StoE,  EPlr), 

where 


S[l,1]  is  the  X  coordinate  and  S[l,2]  is  the  Y  coordinate  of  point  P|  of 

S, 

V[l.1]  is  the  X  coordinate,  V[l,2]  is  the  Y  coordinate  of  the  1th  Voronoi 
point  and  V[l,3]  is  a  one  bit  flag  that  distinguishes  the  Voronoi  points 
of  the  nearest  point  diagram  from  those  of  the  farthest  point  diagram, 

E[l,1]  and  E[l,2]  point  to  the  two  Voronoi  points  of  V  that  determine 
edge  I.  Note  that  for  a  ray  one  point  will  be  from  the  nearest  point 
diagram  and  one  will  be  from  the  farthest  point  diagram.  Furthermore, 
we  store  the  edges  sorted  counterclockwise  to  make  it  convenient  to 
determine  the  Voronoi  polygons  associated  with  each  point  J  of  S.  The 
first  few  edges  of  E  define  the  edges  of  the  nearest  (and  farthest) 
Voronoi  polygons  of  point  1  of  S,  the  succeeding  edges  of  E  define  the 
Voronoi  polygon(s)  for  point  2  of  S,  anti  so  on. 

EPtr[l,1]  points  to  the  first  edge  in  array  E  of  the  nearest  point  Voronoi 
polygon  for  point  I  of  S.  Similarly,  EPtr[l,2]  points  to  the  first  edge  in 
array  E  of  the  farthest  point  Voronoi  polygon  for  point  I  of  S  (if  it 
exists). 

VtoS[l,1],  VtoS[l,2],  and  VtoS[l,3]  point  to  the  three  points  of  S  that 
determine  Voronoi  point  I. 

The  representation  for  the  dual  nearest  and  farthest  point  diagrams  is  equivalent  to 
the  representation  for  the  nearest  and  farthest  point  Voronoi  diagrams. 


4.1.3.  Planar  Voronoi  Diagram  Algorithm 

We  here  combine  the  tools  of  the  previous  sections  to  produce  an  0(N  log  N)  time 
algorithm  for  constructing  a  Voronoi  diagram  of  a  set  S  of  N  planar  points.  The 
algorithm  takes  advantage  of  the  fact  that  the  Voronoi  points  of  the  nearest  point 
diagram  can  be  represented  by  a  set  of  circles  that  each  (i)  pass  through  three  of 
the  N  points  of  S,  and  (ii)  do  not  contain  any  of  these  N  points  in  the  interior.  As 
suggested  in  Section  3.4,  when  we  have  a  problem  that  is  expressed  in  terms  of 
circles,  it  may  be  profitable  to  try  to  apply  inversion  to  obtain  an  equivalent  problem 
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that  involves  linear  components  rather  than  circular  components.  In  this  case  (as 
for  the  union  of  clinks  in  Section  3.2)  wo  must  first  embed  the  planar  problem  into 
three  dimensions  and  express  the  new  problem  in  terms  of  spheres  before  we  can 
profitably  apply  inversion.  This  is  because  inversion  transforms  circular  (or 
spherical)  components  into  linear  components  only  if  the  circle  (or  sphere)  passes 
through  the  center  of  inversion.  We  need  the  extra  degree  of  freedom  that  we 
obtain  by  embedding  in  a  higher  dimension  to  satisfy  this  condition.  The  linear 
(component)  problem  that  we  obtain  turns  cut  to  be  simple:  construct  the  convex 
hull  of  N  transformed  points.  Furthermore,  we  can  also  obtain  the  farthest  point 
diagram  from  the  same  convex  hull. 

Algorithm  for  Construction  of  a  Planer  Voronoi  Diagram 

l.lel  S  be  a  set  of  N  planar  points  located  in  the  xy  plane  of 

three-space.  Pick  a  point  P  in  three-space  that  is  not  in  the  xy  plane 
15 

• 

2.  Choose  any  radius  of  inversion  R  >  0  and  then  invert  the  N  points  of  S 
with  respect  to  point  P  and  radius  R.  (Section  3.2.3  describes 
inversion.)  Call  this  new  set  of  N  points  S'. 

3.  Construct  the  convex  hull  of  the  points  in  S'  in  0(N  log  N)  time  (by  the 
^algorithm  of  Preparata  and  llong  [S3]).  All  N  of  tire  points  of  S’  will  be 

on  (ho  convex  hull  because  inversion  about  P  maps  all  points  of  the  xy 
plane  to  a  sphere  with  P  at  the  apex.  (Sec  Figure  4-5.)  Let  f  be  the 
number  of  faces  on  the  convex  hull  and  let  the  edge  (if  any)  joining 
faces  Fj  and  Fj  be  denoted  Ejj. 

4.  Each  of  tlie  f  faces  Fj  of  the  convex  hull  determines  a  plane  in 
three-space.  Invert  these  f  planes  (with  respect  to  center  of 


Although,  mathematically,  any  point  P  outbid*  the  xy  plane  will  work,  we  may  encounter  excessive  round-off 
error  in  a  computer  if  P  is  badly  chosen.  VVc  warn  to  choose  x  and  y  coordinates  that  (approximately)  center  P 
over  the  convex  hull  of  Ihe  N  points  and  choose  ihc  a  coordinate  so  that  it  is  not  too  close  to  the  xy  plane  (which 
duster*  the  transformed  point*  around  P)  and  also  not  too  far  away  from  the  xy  plane  (which  make*  all  of  the 


transformed  points  approximately  eoplanar).  If  r  , 
coordinates  among  all  of  the  N  planar  points  then 
Pz  ■  mo*(  (hntax  .  ,nW|X  ,min, 

we  can  choose  a  good  pomt  P  m  0(N)  time. 


*mir 

let 


'max’ 
,  *  (*, 


ymtn  *• 

max  + 


■  max- 


the  max  and  mm  x  and  y 
py  *  (ymax  +  ymm^2*  and 
min'  *iwix  "*  *;r«n  in  °<N> 
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inversion  P  and  radius  P.)  to  obtain  f  spheres  that  intersect  the  xy 
plane  in  f  circles.  The  centers  of  these  circles  are  the  Voronoi  points 
Vj.  To  distinguish  nearest  and  farthest  Voronoi  points  we  perform  the 
following  simple  test: 

The  plane  of  face  Fj  determines  two  half-spaces,  one 
that  contains  the  entire  convex  hull  and  one  that  contains 
none  of  it.  Let  half-space  llj  be  the  one  that  contains  the 
convex  hull.  If  llj  contains  point  P  then  Vj  is  a  Voronoi  point 
of  the  nearest  point  Voronoi  diagram.  Otherwise,  Vj  is  a 
Voronoi  point  of  the  farthest  point  diagram. 

5.  We  obtain  the  Voronoi  points  Vj  from  the  faces  Fj  of  the  convex  hull, 
but  to  construct  the  remainder  of  the  Voronoi  diagram  we  must 
examine  the  edges.  Each  edge  Ejj  of  the  convex  hull  corresponds  to  a 
segment  of  the  nearest  point  diagram,  a  segment  of  the  farthest  point 
diagram,  or  a  ray  (for  both  diagrams).  To  determine  for  an  edge  Ejj 
which  of  these  three  possibilities  is  true  we  use  the  following  rules: 

-  If  V;  and  Vj  are  both  nearest  Voronoi  points,  then  there  is  a  line 
segment  connecting  Vj  and  Vj  in  the  nearest  point  diagram. 

-  If  Vj  and  Vj  are  both  farthest  Voronoi  points,  then  there  is  a  line 
segment  connecting  Vj  and  Vj  in  the  farthest  point  diagram. 

-  If  Vj  is  a  closest  Voronoi  point  and  Vj  is  a  farthest  Voronoi  point, 
then  Vj  and  Vj  determine  a  ray  in  both  the  nearest  and  farthest 
point  Voronoi  diagrams.  The  points  Vj  and  Vj  determine  a  line, 
and  the  desired  ray  for  the  nearest  point  Voronoi  diagram  is  the 
part  of  that  line  that  starts  at  point  Vj  and  does  not  include 
point  Vj.  The  ray  starting  at  point  Vj  that  does  not  include  point 
Vj  is  for  the  farthest  point  Voronoi  diagram. 

Although  it  is  clear  that  the  above  algorithm  requires  only  0(N  log  N)  time  and  O(N) 
storage,  it  is  not  immediately  obvious  that  it  actually  constructs  the  nearest  (or 
farthest)  point  Voronoi  diagram.  We  must  explain  (1)  why  the  centers  of  the  circles 
(generated  in  Stop  4  above)  arc  the  Voronoi  points  and  (2)  why  the  connection 
rules  (Step  5)  for  Voronoi  points  work.  We  wilt  no w  show  this  for  the  case  of  the 
nearest  point  diagram.  The  argument  for  the  farthest  point  diagram  is  similar. 

Theorem  20:  The  centers  of  the  circles  generated  in  Step  4  of  the 
above  algorithm  ore  the  Voronoi  points. 
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Figure  4-5:  Planar  Voronoi  diagram  and  corresponding  convex  hull. 

Proof:  The  proof  is  in  two  parts:  ( 1 )  all  of  the  points  generated  in  Step 
4  are  Voronoi  points,  and  (2)  all  of  the  Voronoi  points  are  generated  in 
Step  4.  To  prove  that  the  circles  generated  in  Step  4  are  centered  at 
the  Voronoi  points  we  must  show  that  (a)  these  circles  each  pass  through 
three  of  the  N  points  of  S  and  (b)  do  not  contain  any  of  the  other  N  -  3 
points  in  the  interior  (Section  4.1.1).  Part  (a)  follows  from  the  fact  that 
inversion  is  involutory  (Section  3.2.3).  We  prove  part  (b)  by 
contradiction.  Assume  that  the  circle  passing  through  points  A,  B,  and  C 
of  S  contains  another  point  Q  f  S  in  its  interior.  This  places  point  Q  inside 
the  sphere  determined  by  points  A,  B,  C,  and  P.  When  we  invert  about 
point  P,  the  point  Q'  is  separated  from  point  P  by  the  plane  determined  by 
points  A',  B',  and  C’.  Since  the  plane  A'B'C'  does  not  determine  a 
half-space  that  contains  point  P  and  also  contains  all  of  the  other  N  -  3 
points  of  S'  it  cannot  be  a  face  of  the  convex  hull  that  determines  a 
Voronoi  point  of  the  nearest  point  Voronoi  diagram. 

To  prove  that  Step  4  generates  all  of  the  Voronoi  points  we  simply  use 
the  reverse  argument.  If  Vj  is  a  Voronoi  point  then  the  circle  for  Vj 
transforms  (by  inversion)  to  a  face  of  the  convex  hull  of  S’.  Let  A,  B,  and 
C  be  the  points  of  S  that  determine  Voronoi  point  V.  The  points  A',  B*,  and 
C*  of  S'  determine  a  plane  that  contains  all  of  the  points  of  S'  because  all 
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N  -  3  other  points  of  S  lie  outside  the  circle  determined  by  points  A,  B, 

and  C.  □ 

Theorem  PU  Step  5  of  the  algorithm  correctly  obtains  the  edges  of  the 
Voronoi  diagram. 

Proof:  The  proof  is  in  two  parts:  (1)  all  of  the  edges  generated  by  Step 
5  are  edges  of  the  diagram  and  (2)  all  edges  of  the  diagram  are 
generated  by  Step  5.  An  edge  Ejj  of  the  convex  hull  that  separates 
(nearest  point)  faces  Fj  and  Fj  maps  to  a  line  segment  between  Voronoi 
points  Vj  and  Vj.  Blit  the  circles  corresponding  to  Voronoi  points  V,  and  Vj 
meet  at  two  of  the  N  points  of  S  because  the  corresponding  faces  Fj  and 
Fj  share  an  edge  Fjj.  This  is  exactly  the  characterization  given  for  edges 
of  the  Voronoi  diagram  in  Section  4.1.1.  Similarly,  the  rays  are  determined 
by  edges  Ejj  where  Vj  is  a  nearest  Voronoi  point  and  Vj  is  a  farthest 
Voronoi  point  (or  vice-versa).  In  this  case,  too,  the  circles  corresponding 
to  Vj  and  Vj  meet  at  two  of  the  N  points  of  S.  □ 

Since  the  Voronoi  points  and  edges  (and  rays)  connecting  the  Voronoi  points  are 

correctly  generated  by  the  above  algorithm,  we  have  just  proven 

Theorem  22:  The  algorithm  constructs  the  Voronoi  diagram  in  0(N  log  N) 
time. 


4.1.4.  Fast  Expectcd-Timc  Algorithms 

The  most  expensive  part  of  the  algorithm  for  construction  of  a  Voronoi  diagram  is 
the  construction  of  the  convex  hull.  If  the  convex  hull  can  be  constructed  in  fast 
expected  time,  then  the  Voronoi  diagram  can  be  constructed  in  fast  expected-time. 
The  O(N)  expected-time  algorithms  of  Bentley  and  Shamos  [16],  Eddy  [39],  or  Floyd 
[40]  do  not  apply  because  their  results  depend  on  a  sublinear  expected  number  of 
points  on  the  convex  hull,  and  for  the  Voronoi  diagram  algorithm  there  are  always  N 
vertices  on  the  convex  hull. 

Bentley,  Wcidc,  and  Yao  [IS],  on  the  other  hand,  describe  how  a  planar  Voronoi 
diagram  can  be  constructed  in  linear  expected-time.  The  only  condition  is  that  the 
probability  density  of  the  underlying  distribution  must  be  bounded  above  and  below 
by  (nonzero)  constants.  The  algorithm  does  not  make  use  of  inversion.  Instead,  it 
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applies  an  extension  of  Weide's  [99]  technique  for  an  0(N)  expected-time  sort  to 
the  planar  Voronoi  diagram  problem. 

4.1.5.  Higher  Dimensions 

The  K-dimensional  Voronoi  diagram  algorithm  is  an  extension  of  the  planar 
algorithm.  We  first  embed  the  M  K-dimensional  points  of  S  in  K+1  -space  and  then 
invert  them  to  N  K+1  -dimensional  points  S'.  We  then  construct  the  convex  hull  of  S' 
and  obtain  the  Voronoi  diagram  by  transforming  the  parts  of  the  convex  hull  back  to 
K-space.  To  transform  back  to  K-space  we  first  invert  each  hyperface  of  the 
convex  hull  to  obtain  a  set  of  K+1  -spheres  whose  intersection  with  K-space  is  a  set 
of  K-spherns.  These  K-spheres  each  pass  through  K+1  points  of  S  and  are 
centered  at  the  Voronoi  points.  We  obtain  the  other  components  of  the 
K-dimensional  Voronoi  diagram  by  connection  rules  similar  to  those  in  Step  5  of  the 
algorithm  in  Section  4.1.3.  For  example,  if  the  K-sphere  for  Voronoi  point  Vj  passes 
through  K  of  the  K+1  points  determining  the  K-sphere  for  Voronoi  point  Vj,  then  we 
draw  a  one-dimensional  edge  between  V;  and  Vj.  If  the  spheres  for  a  set  of  three  or 
more  Voronoi  points  share  K-1  points  of  S,  then  we  draw  a  two-dimensional  edge 
between  the  Voronoi  points  of  that  set.  (A  two-dimensional  edge  between  L  points 
is  a  convex  polygon  with  L  vertices.)  The  rules  for  three  and  higher  dimensional 
edges  are  similar.  The  time  complexity  of  the  K-dimensional  Voronoi  diagram 
algorithm  is  dominated  hy  the  time  to  construct  a  K+1 -dimensional  convex  hull  of  N 
points.  (See  Section  1.1.1  for  references  to  several  convex  hull  algorithms.) 

4.2.  Spherical  Nearest  and  Farthest  Point  Voronoi  Diagrams 

Voronoi  diagrams  are  useful  for  solving  several  closest  or  farthest  point 
geographic  problems.  If,  hov/ever,  the  area  covered  by  the  points  is  large,  then  we 
must  take  the  curvature  of  the  earth  into  consideration.  The  most  obvious 
approximation  to  use  for  the  earth  is  a  sphere.  Nearest  and  farthest  point  Voronoi 
diagrams  on  a  sphere  are  defined  in  a  manner  analogous  to  their  planar  counterparts 
and  the  algorithms  for  constructing  them  provide  an  interesting  comparison  with 
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those  for  the  planar  case.  For  example,  we  obtain  the  farthest  point  Voronoi 
diagram  of  a  set  of  N  spherical  points  S  by  simply  applying  a  nearest  point  algorithm 
to  a  set  S'  of  N  points  diametrically  opposite  the  points  of  S.  There  are  two 
different  methods  for  constructing  these  diagrams  in  0(N  log  N)  time.  One  involves 
an  intersection  of  half-soaccs  and  the  other  obtains  the  dual  of  the  Voronoi  diagram 
from  the  convex  hull  of  the  spherical  points.  We  will  describe  only  the  second 
algorithm  because  it  is  simpler. 

We  take  advantage  of  the  fact  that  the  spherical  Voronoi  diagram,  as  well  as  the 
planar  Voronoi  diagram,  can  be  expressed  in  terms  of  a  set  of  circles:  the  Voronoi 
points  are  the  centers  of  the  circles  (on  the  sphere)  that  (i)  pass  through  three  of 
the  N  spherical  points,  and  (ii)  do  not  contain  any  of  the  other  N  -  3  spherical  points. 
As  before,  pairs  of  circles  that  siiare  two  of  the  N  points  determine  the  edges  of  the 
diagram.  We  would  like  to  express  this  problem  in  terms  of  linear  components  rather 
than  circular  components. 

One  approach  is  to  construct  a  spherical  analog  of  the  formula  for  the  planar 
Voronoi  diagram  algorithm:  embed  to  a  four-dimensional  sphere,  apply  (spherical) 
inversion  (with  respect  to  a  suitable  point  P  of  the  four-sphere),  and  construct  the 
(spherical)  convex  hull  of  the  transformed  points.  Although  this  approach  can 
actually  be  made  to  work,  it  does  not  give  us  a  problem  that  involves  linear 
components.  A  belter  approach  is  to  embed  the  spherical  Voronoi  diagram  problem 
(which  is  a  spherical  two-space  problem)  into  Euclidean  three-space.  The  circles 
that  define  the  spherical  Voronoi  diagram  determine  the  planes  that  bound  the  faces 
of  the  (Euclidean  three-dimensional)  convex  hull  of  the  N  points.  This  convex  hull 
can  be  converted  readily  into  the  dual  of  the  spherical  Voronoi  diagram  —  the 
spherical  Delaunay  diagram  --  and  can  be  constructed  ii  0(N  log  N)  time.  Given  the 
dual,  the  Voronoi  diagram  can  be  produced  in  only  0(N)  additional  time. 
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Algorithm  for  Spherical  Vorouoi  Diagram 

1.  Let  S  be  a  set  of  N  >  3  points  on  the  surface  of  a  sphere  such  that  no 
four  points  are  co-circular.  Construct  the  convex  hull  of  the  points  of 
S  (treating  them  as  N  points  in  Euclidean  three-space)  in  0(N  log  N) 
time  by  the  algorithm  of  Prcparnta  and  Hong  [83]. 

2.  For  each  face  Fj  of  the  convex  hull  there  is  a  corresponding  Voronoi 
point  Vj  on  the  surface  of  the  sphere  that  is  equidistant  from  the 
vertices  of  face  Fj.  (Actually  there  are  two  such  points  —  V;  and  the 
point  diametrically  opposite  Vj.  Choose  the  point  that  is  closest  to  the 
vertices  of  face  Fj.) 

3.  For  each  pair  of  faces  Fj  and  Fj  that  share  on  edge  Ejj  construct  an 
arc  of  a  great  circle  that  connects  points  Vj  and  Vj.  Since  there  ore 
only  O(N)  faces  and  edges  in  the  convex  hull  we  can  do  this  in  O(N) 
time. 

The  diagram  that  the  above  algorithm  constructs  is  the  spherical  nearest  point 
Voronoi  diagram  because  each  Voronoi  point  Vj  corresponding  to  face  Fj  is  not  only 
equidistant  from  the  vertices  of  face  Fj  but  is  also  closer  to  these  three  points  than 
the  other  N-3  points  of  set  S.  As  in  the  case  of  the  planar  Voronoi  diagram,  we  have 
a  set  of  circles  that  ensure  the  properties  of  the  Voronoi  diagram. 


4.3.  Nearest  and  Farthest  Edge  Diagrams 

Drysdale  and  I  ne  f371  describe  the  construction  of  Voronoi  dianrams  of  line 
segments  (and  other  geometrical  objects)  in  0(N  c^°3  N)l/‘)  tjme  (which  they  later 
improved  to  0(N  (log  N)^)  time)  and  Kirkpatrick  [59]  has  reduced  this  time  to 
0(N  log  N).  A  special  case  of  this  problem  is  the  nearest  (respectively  farthest) 
edge  diagram  for  a  convex  N-gon.  This  diagram  is  a  tcsselation  of  the  plane  into  N 
polygonal  regions  such  that  each  region  i  is  the  set  of  all  points  nearest  to 
(respectively  farthest  from)  edge  i  of  the  N-gon.  Figure  4-6  illustrates  a  nearest 
edge  diagram. 


The  problem  of  constructing  a  nearest  edge  diagram  is  presented  by  Shomos 
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[89]  as  problem  P0L9  in  bin  workbook.  One  interesting  application  is  that  once  we 
have  constructed  the  nearest  edge  diagram,  we  can  solve  the  problem  of 
constructing  the  greatest  incircle  cf  a  convex  polygon  (problem  P0L1G)  in  0(N)  time. 
Both  Preparata  and  Lee  describe  0(N  log  N)  time  solutions  to  this  problem  that  are 
not  based  on  geometric  transforms  [77,  79].  (They  call  this  problem  "Medial  axis  of 
a  convex  polygon.")  In  this  section  we  describe  an  0(N  log  N)  time  solution  that  is 
based  on  the  simple  geometric  transforms  of  (1)  embedding  in  a  higher  dimension 
and  (2)  orthographic  projection. 

One  of  the  main  differences  between  the  nearest-edge  diagram  and  the  others 
described  in  this  chapter  is  that  the  elements  defining  the  diagram  are  edges  rather 
than  points.  We  can,  however,  still  express  it  in  terms  of  circles;  since  each  vertex 
V  of  the  diagram  is  equidistant  from  the  three  nearest  edges  of  the  N-gon,  V  is  the 
center  of  a  circle  that  is  tangent  to  three  edges  (but  does  not  intersect  or  contain 
any  of  the  other  N  -  8  edges).  This  suggests  that  we  might  try  to  apply  the  formula 
that  we  used  Tor  the  Euclidean  planar  Voronoi  diagram  of  N  points;  embed  in 
Euclidean  three-space,  invert  the  edges  with  respect  to  a  point  P  that  is  not  in  the 
xy  plane  (producing  a  connected  set  of  circular  arcs),  and  construct  the  convex  hull 
of  the  transformed  elements.  This  would  work  well  if  we  had  a  fast  algorithm  for 
constructing  the  convex  hull  of  a  set  of  (connected)  circular  arcs  in  three-space.1® 
There  is,  however,  another  way  to  embed  this  problem  in  three-space  that  produces 
a  three-dimensional  problem  that  involves  only  linear  components. 

The  N  circular  arcs  all  determine  circles  that  pass  through  point  P.  For  any  one  of 
these  N  circles  Cj  there  is  an  infinite  number  of  spheres  that  pass  through  (all  of  the 
points  of)  C,.  We  can  thus  represent  the  convex  N-gon  (in  the  xy  plane)  as  a  set  of 
N  spheres  (in  three-space)  each  of  which  passes  through  point  P  and  still  has  one 
degree  of  freedom.  Ilow  should  the  spheres  he  chosen? 


* 
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Figure  4-3:  Nearest  edge  diagram  of  a  convex  polygon. 

Since  inversion  maps  spheres  to  planes  and  the  interior  of  the  spheres  to 
heif-spaccs.  we  con  represent  the  nearest  edge  diagram  by  an  intersection  of 
half-spaces.  The  problem  of  choosing  the  spheres  now  becomes  a  problem  of 
choosing  the  half-spaces.  For  each  haif-spaee  the  degree  of  freedom  is  the  angle 
that  its  boundary  makes  with  the  xy  plane.  Since  the  edges  of  the  nearest  edge 
diagram  lie  on  the  angular  bisectors  of  adjacent  sides  of  the  N-gon  we  choosa 
half-spaces  whose  boundaries  meed  directly  above  these  angular  bisectors.  This 
gives  us  the  following  algorithm: 
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Algorithm  for  Nearest  (Respectively  Farthest)  Edge  Diagram  for  a  Convex 
Polygon 

Input:  N  --  number  of  vertices  in  the  convex  polygon.  X[1:N],  Y[1:N]  —  x  and  y 
coordinates  of  vertices  of  the  M-gon  (counterclockwise  order). 

Output:  Nearest  Edge  Diagram  The  representation  is  similar  to  the  representation  for 
an  intersection  of  half-spaces  (Section  3.1). 

Time:  0(N  log  N),  Space:  0(N). 

1 .  Let  the  convex  polygon  lie  in  the  xy  plane  of  3-space.  For  each  edge 
ej  of  the  polygon  construct  the  unique  plane  Pj  that  (a)  contains  that 
edge,  (b)  makes  a  45  degree  angle  with  the  xy  plane,  and  (c)  lies 
above  the  polygon  (rather  than  below  it).  For  a  nearest  edge  diagram, 
let  h|  be  the  half-space  that  lies  below  plane  p;.  For  a  farthest  edge 
diagram,  let  hj  lie  above  plane  Pj. 

2.  Intersect  the  N  half-spaces  hj  in  0(N  log  N)  time  (Section  3.1 ). 

3.  Project  the  intersection  to  the  xy  plane.  (This  amounts  to  throwing  out 
the  z  coordinates  of  the  vertices  of  the  intersection.) 

The  above  algorithm  does  not  make  any  distance  measurements  to  construct  the 
nearest  (or  farthest)  edge  diagram.  Instead,  it  relies  on  the  symmetry  induced  in 
Step  1  by  constructing  all  planes  Pj  at  the  same  angle  (45  degrees)  from  the  xy 
plane.  We  can  generalize  this  to  weighted  distances  Wj  from  the  edges  by  simply 
letting  the  slopes  of  the  planes  Pj  be  set  to  the  weights  Wj. 

4.4.  Summary 

We  have  described  three  types  of  diagrams,  planar  and  spherical  Voronoi 
diagrams  of  sets  of  points,  and  (nearest  and  farthest)  edge  diagrams  for  a  convex 
polygon.  Since  each  of  these  problems  involved  Euclidean  distance  between  the 
elements  defining  the  diagram,  we  found  it  useful  to  express  each  problem  in  terms 
of  circles.  We  then  embedded  the  problem  in  a  higher  dimension  and,  when 
necessary,  expressed  it  in  terms  of  (carefully  chosen)  spheres  and  applied 
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inversion  to  obtain  an  equivalent  problem  expressed  in  terms  of  linear  components.. 

For  the  planar  Vorouoi  diagram  problem  wc  embedded  the  plane  into 
three-dimensions  and  applied  inversion  to  obtain  a  convex  hull  problem.  The  circles 
defining,  the  planar  Vorouoi  diagram  are  centered  at  the  Voronoi  points  and  each 
pass  through  three  of  the  N  points  but  do  not  contain  any  of  the  other  N  -  3  points  in 
their  interiors.  These  circles  become  spheres  when  embedded  into  three  dimensions 
and,  when  inverted,  become  the  planes  that  bound  the  faces  of  the  convex  hull. 

We  were  able  to  solve  the  spherical  Voronoi  diagram  problem  directly  as  a  convex 
hull  problem.  This  is  because  the  circles  defining  the  diagram  pass  through  three  of 
the  N  spherical  points  but  do  not  contain  any  of  the'  other  N  -  3  points.  When  the 
spherical  Voronoi  diagram  problem  (which  is  a  spherical  two-space  problem)  is 
embedded  into  Fuc.lidoan  three-space,  these  circles  become  the  planes  that  bound 
the  convex  hull  of  the  N  points. 

We  constructed  the  nearest  (respectively  farthest)  edge  diagram  of  a  convex 
polygon  by  embedding  into  three-space  and  intersecting  half-spaces.  Each  of  the 
circles  that  define  the  nearest  edge  diagram  is  tangent  to  three  of  the  sides  of  the 
convex  N-gon  and  does  not  contain  any  part  of  the  other  N  -  3  sides  in  its  interior. 
By  embedding  into  Euclidean  three-space,  applying  inversion,  redefining  the  problem 
in  terms  of  spheres,  and  then  applying  inversion  again  we  obtain  a  problem  of 
intersecting  half-spaces.  But,  as  we  saw  in  Section  3.1,  the  intersection  of 
half-spaces  is  solved  by  constructing  the  convex  hull  of  a  set  of  points. 

.  In  the  next  chapter  we  will  find  evert  more  uses  for  convex  hulls. 
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5.  Searching  Tesselations 

In  this  chapter  we  demonstrate  how  a  search  of  a  tcssclation  (Section  1.1.4) 
arises  in  both  linear  programming  and  computing  the  diameter  of  a  set  of  points. 
Both  of  these  problems  invite  the  use  of  an  orthographic  projection  to  reduce  a 
K-dlmensionnl  problem  to  a  K-1 -dimensional  problem  and  furthermore  provide 
interesting  applications  of  the  point  /  flat  duality  transform  (Section  3.1 .3.3).  The 
diameter  of  a  set  of  points  also  provides  another  application  of  the  convex  hull  of  a 
set  of  points. 

5.1.  Linear  Programming 

Linear  programming  is  an  important  technique  for  optimizing  a  linear  function 
subject  to  a  set  of  several  linear  constraints,  if  there  are  K  variables  and  N 
constraints,  we  may  interpret  each  constraint  as  a  half-space  in  K-space  and  the 
feasible  region  satisfying  all  of  the  constraints  as  the  polytope  that  is  the 
intersection  of  N  K-dimcnsional  half-spaces.  One  of  the  vertices  of  the  resulting 
polytope  is  an  optimal  solution  for  the  linear  program.  The  linear  programming 
problem  is  to  find  this  vertex  as  quickly  as  possible. 

The  standard  method  of  solving  linear  programming  problems  is  the  simplex  method 
(and  its  variants)  [28],  In  the  worst  case,  however,  the  simplex  method  will  require 
exponential  time  [60,  55].  Kelly  [OS]  has  shown  that  for  a  model  of  linear 
programming  with  M  relevant  constraints  in  two  variables  chosen  independently  from 
a  given  distribution,  the  expected  number  of  iterations  is  0(N).  Since  each  iteration 
costs,  in  this  case,  0(N)  time,  the  expected  time  for  linear  programming  in  his 
two-dimensional  model  is  O(N^).  In  general,  however,  the  expected-time  of  the 
simplex  method  has  not  been  adequately  analyzed,  although  empirical  results 
indicate  that  it  may  be  bounded  above  by  a  low  degree  polynomial  in  K  and  N  [33]. 

There  are  alternatives  to  the  simplex  method.  Khachian  [44]  has  produced  an 
algorithm  for  solving  linear  programming  with  integral  coefficients  that  costs  only 
polynomial  time  (in  the  size  of  the  input)  in  the  worst  case.  His  method,  however, 
depends  strongly  on  the  fact  that  the  coefficients  are  integers.  Another  approach 
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is  to  intersect  the  N  K-tiimensional  half-spaces  to  construct  explicitly  the  feasible 
region  and  then  evaluate  the  objective  function  at  each  vertex.  For  K  =  2  or  K  *  3 
variables  we  may  construct  the  intersection  and  solve  the  linear  programming 
problem  in  0(N  log  N)  worst  case  time  and  (wh.cn  the  expected  number  of 
nonredumlont  constraints  is  0(NP)  for  some  p  <  1)  0(N)  expected-time 
[94,  106,  84]  This  is,  in  the  worst  case,  better  than  the  simplex  method,  which  may 
take  O(N^)  time  in  these  cases.  It  is  not,  however,  generally  considered  advisable 
to  construct  explicitly  the  entire  feasible  region  when  we  need  only  the  vertex 
corresponding  to  the  solution.  The  simplex  method  certainly  avoids  that  problem. 
We  now  describe  an  approach  toward  a  better  solution. 


Figure  5-1:  Transform  of  LP  problem  to  vertical  line  and  convex  hull. 

Dantzig  [28]  describes  an  alternate  interpretation  of  the  simplex  method  that  is 
based  on  the  point  /  flat  duality  transform.  The  polytope  obtained  by  intersecting 
the  N  K-dimensional  half-spaces  transforms  to  the  convex  hull  of  N  points  in 
K-space  and  the  (linear)  objective  function  transforms  to  a  vertical  line  (in  the  K 
coordinate)  in  K-space.  (See  Figure  5-1.)  The  linear  programming  problem  itself  is 
transformed  to  the  problem  of  determining  which  face  of  the  convex  luiil  this  vertical 
line  intersects.  When  we  orlhogrophically  project  the  convex  hull  to  a  tesselation 
and  the  verlical  line  to  a  point  in  K-1 -space,  wc  obtain  a  problem  of  locating  a  point 
in  a  tesselation.  (See  Figure  5-2  for  an  illustration  of  the  two-dimensional  case.) 
This  docs  not  immediately  lead  to  a  (provably)  faster  linear  programming  algorithm 
than  the  simplex  method,  but  it  does  provide  another  way  to  approach  the  linear 
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programmiiuj  problem. 


Figure  5-2:  Transform  of  vertical  line  and  convex  hull  to  a  point  and  a  tesselation. 


5.2.  Diameter  of  a  Set  of  Points 

The  diameter  of  n  set  of  points  is  the  distance  between  the  two  farthest  points. 
This  quantity  often  arises  in  problems  of  cluster  analysis  [51]  because  a  set  of 
points  that  are  all  near  each  other  makes  a  better  cluster  than  a  set  of  points  that 
are  spread  far  apart.  The  straightforward  way  to  determine  the  diameter  of  N 
points  is  to  compute  all  0(N"-)  interpoint  distances  and  return  the  maximum. 
Depending  on  the  metric  and  the  number  of  dimensions,  however,  there  may  be  much 
faster  algorithms. 

For  K  *  1  dimension  oil  N  points  lie  on  a  line  and  the  diameter  is  the  distance 
between  the  points  with  maximum  and  minimum  coordinate,  which  can  easily  be  found 
in  O(N)  time.  <ln  fact,  F3N/2  *  2l  comparisons  are  necessary  and  sufficient 
[76]).  For  two  or  more  dimensions  the  choice  of  metric  is  important.  The  L-j  or 
diameter  of  N  points  in  K  dimensions  con  be  easily  computed  in  0(2KN)  or  0(KN)  time, 
respectively,1^  but  the  Euclidean  case  may  be  more  difficult.  For  two  dimensions 


17Tho  diameter  of  a  set  is  determined  bv  Stic  most  extreme  oointj  in  each  of  the  direction*  determined  bv  the 

'  ‘  • 

feces  of  the  unit  "inhere".  Since  the  unit  inhere  for  the  L.  metric  hai  2  faces,  and  the  unit  sohere  for  the  toe 
metric  has  2K  faces,  we  can  find  the  extreme  points  in  0(2*M)  and  O(KiN)  time,  respectively. 
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the  best  known  Euclidean  diameter  algorithms  run  in  worst-case  time  0(N  log  N) 
[89,  91].  For  three  dimensions  Yao  [10?]  tins  produced  an  0(N1-®)  time  algorithm 
and  for  K  dimensions  0(N^  *  time,  where  oT-CK)  =  2"^+0. 


We  will  present  an  0((N  +  K)  log  N)  time  algorithm  for  the  three-dimensional  case, 
where  K  is  the  number  of  pairs  of  antipodal  vertices  on  the  convex  hull.  To  achieve 
this  time  we  apply  a  point  /  flat  duality  combined  with  orthographic  projection  to  the 
components  of  the  convex  hull  of  the  set  of  N  points  to  obtain  a  problem  of  locating 
points  in  an  outerplanar  straight-line  graph.  In  Appendix  II  we  describe  a 
relationship  between  the  Euclidean  diameter  and  an  empty-intersection  problem  that 
may  lead  to  an  ft(N  log  N)  time  lower  bound  for  the  diameter  problem.  In  the 
following  sections  we  present  0(N  log  N)  time  two-dimensional  Euclidean  diameter 
algorithms,  the  0((N  +  K)  log  N)  time  three-dimensional  algorithm,  and  then  discuss 
fast  expected-time  algorithms,  approximation  algorithms,  higher  dimensions, 
applications,  and  some  unsolved  problems. 

5.2.1.  Diameter  in  Two  Dimensions 

Shamos  [91]  describes  an  0(N  log  N)  time  algorithm  for  computing  the  diameter  of 
N  points  in  the  plane.  The  algorithm  that  wc  present  is  essentially  equivalent  to  his, 
but  it  is  expressed  so  that  it  generalizes  to  a  fast  three-dimensional  algorithm.  We 
first  present  a  theorem  that  reduces  our  search  for  the  diametrical  pair  of  points  to 
the  convex  hull. 

Theorem  23:  (Hocking  and  Young  [52],  p.  207)  The  Euclidean  diameter 
of  a  set  of  points  S  is  determ. ned  by  two  points  on  the  convex  hull  of  S. 

If  all  N  of  the  points  of  S  are  on  the  convex  hull,  then  we  have  not  reduced  the 
size  of  the  problem.  We  have,  however,  simplified  it  by  reducing  the  problem  of 
computing  the  diameter  of  a  set  of  points  to  the  problem  of  computing  the  diameter 
of  a  convex  polygon.  For  our  new  problem  we  have  the  following  theorem: 

Theorem  24:  (Yaglom  and  Bollyanskii  [100],  p.  9)  The  diameter  of  a 
convex  figure  is  the  maximum  distance  between  parallel  lines  of  support 
of  this  figure. 


Figure  5-3:  Convex  hull  and  parallel  lines  of  support. 

In  Figure  5-3  we  illustrate  two  (parallel)  lines  of  support  and  L2.  In  general,  a 
line  of  support  passes  through  (at  least)  one  boundary  point  of  a  figure  and  lies 
entirely  on  one  side  of  that  figure.  Pairs  of  points  (of  the  figure)  on  opposite 
parallel  lines  of  support  are  called  antipodal  points.  In  Figure  5-3  points  A  and  E  and 
points  B  and  E  are  antipodal.  We  are  interested  in  the  antipodal  pairs  of  vertices 
determined  hy  the  lines  of  support  for  the  convex  hull  of  S  because  one  of  these 
pairs  determines  the  diameter  of  S.  Our  next  theorem  gives  a  bound  on  the  number 
of  pairs  that  we  will  have  to  examine. 

Theorem  25:  For  a  convex  polygon  of  N  vertices  there  are  only  0(N) 
antipodal  pairs  of  vertices. 

Proof:  As  we  rotate  parallel  lines  of  support  L  and  M  about  the  convex 
polygon,  the  antipodal  vertices  determined  by  l  and  M  change  only  when 
either  L  or  M  becomes  coincident  with  one  of  the  N  sides.  We  may  thus 
generate  all  of  the  pairs  of  antipodal  pairs  of  vertices  by  recording  all 
antipodal  pairs  when  either  L  or  M  contains  a  side  of  the  polygon.  If  line  L 
contains  a  side,  then  it  passes  through  two  vertices  and  similarly,  its  line 
of  support  M  will  pass  through  at  most  two  vertices.  (In  fact,  only  when 
the  polygon  has  parallel  sides  can  both  L  and  M  simultaneously  pass 
through  two  vertices.)  There  arc  therefore  at  most  four  antipodal  pairs  of 
vertices  generated  each  time  a  line  of  support  passes  through  a  side  of 
the  polygon.  Since  there  are  only  N  sides  of  the  polygon,  there  are  only 
0(N)  antipodal  pairs  of  vertices,  u 


24  December  19^9. 


Geometric  Transforms 


v 

i 


PAGE  86 


We  now  have  enough  information  to  outline  our  two-dimensional  diameter 
algorithm: 

Outline  of  Two-Dimensional  Diameter  Algorithm 

1.  Construct  the  convex  hull  of  the  N  planar  points. 

2.  Generate  the  0(N)  antipodal  pairs  of  vertices  from  the  lines  of  support 
of  the  convex  hull. 

3.  Compare  the  distances  between  each  pair  of  antipodal  vertices  and 
report  the  maximum  as  the  diameter. 

We  can  construct  the  convex  hull  (Step  1)  in  0(N  log  N)  time  [43]  and  easily 
compare,  the  distances  between  antipodal  pairs  of  vertices  (Step  3)  in  only  O(N) 
time.  We  yuis  have  only  to  determine  how  fast  we  can  generate  the  O(N)  pairs  of 
antipodal  vertices  (Step  2).  In  the  remainder  of  this  section  we  show  two 
algorithms  for  generating  them  in  0(N)  time,  making  the  total  time  for  our  diameter 
algorithm  0(N  log  N). 

Shamos  [91]  describes  in  detail  how  to  generate  the  0(N)  pairs  of  antipodal 
vertices  of  a  convex  polygon  in  C(N)  time.  After  finding  the  first  pair  he  generates 
the  other  pairs  in  a  counterclockwise  scan  about  the  polygon,  maintaining  parallel 
lines  of  support  at  all  times.  For  example,  in  Figure  5-3  line  passes  through 
vertices  A  and  B  and  parallel  line  of  support  !_£  passes  through  vertex  E.  We  can 
rotate  counterclockwise  about  vertex  B  and  1 2  about  vertex  E  until  either  L-j 
contains  side  DC  or  L2  contains  side  EF.  We  determine  which  of  the  two  possibilities 
occurs  first  by  comparing  the  slopes  of  sides  EF  and  BC.  In  this  case  line  L2  will 
meet  side  EF  before  line  L-j  meets  CC  because  the  slope  of  EF  is  less  than  the  slope 
of  BC.  This  means  that  vertices  B  and  F  ore  antipodal  and  that  we  will  begin 
rotating  I2  about  vertex  F  rather  than  E.  We  continue  this  procedure  until  the 
parallel  lines  of  support  I}  and  L2  have  traversed  the  entire  convex  polygon  (and 
have  thus  generated  all  of  the  antipodal  pairs  of  vertices). 
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We  can  modify  Slinmos'  algorithm  so  that  it  generalizes  easily  to  the 
three-dimensional  case.  The  important  feature  that  we  extract  from  his  algorithm  is 
that  when  we  perform  the  0(N)  time  scan  around  the  convex  hull,  the  only 
comparisons  that  we  must  make  are  comparisons  of  the  slopes  of  the  sides  of  the 
convex  hull.  That  is,  the  x  and  y  coordinates  of  the  vertices  do  not  matter  since  we 
compare  only  the  slopes  of  the  sides  of  the  polygon.  This  insight  leads  to  a 
one-dimensional  interpretation  of  the  (originally)  two-dimensional  problem. 


Figure  5-4:  Transform  of  a  convex  hull  to  a  line. 

We  illustrate  the  transform  in  Figure  5-4.  The  first  step  is  to  divide  the  convex 
polygon  into  two  parts,  UPPER  and  LOWER.  This  ensures  that  when  one  of  two 
parallel  lines  of  support  meets  the  polygon  at  an  UPPER  side,  the  other  line  of 
support  meets  it  a  LOWER  side.  Wc  may  define  the  transform  of  an  UPPER  side  of 
the  convex  hull  as  the  slope  of  the  UPPER  line  of  support  that  contains  it.  We  may 
also  transform  an  UPPER  vertex  V  to  the  set  of  slopes  of  all  UPPER  lines  of  support 
that  pass  through  V.  The  transform  for  LOWER  sides  and  vertices  is  similar. 
Furthermore,  since  we  consider  the  leftmost  and  rightmost  vertices  (A  and  E  in 
Figure  5-4)  to  lie  in  both  the  UPPER  and  LOWER  sets,  they  each  have  both  an  UPPER 
and  a  LOWER  transform.  We  now  describe  the  transform  algebraically. 

The  transform  maps  the  UPPER  set  of  sides  of  the  convex  hull  to  a  set  of  points 
on  the  line  and  the  LOWER  set  of  sides  to  another  set  of  points  on  the  tine.  The 
mopping  is  simply  a  point  /  flat  duality  followed  by  an  orthographic  projection 
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y  =  nix  +  b  -*  (m,b)  -»  (m)  (14) 

where  "y  =  nix  +  b"  is  the  line  determined  by  a  side  of  the  convex  hull  and  "(m)"  is 
the  one-dimensional  point  to  which  the  side  maps.  The  transform  of  a  vertex  V  of 
the  convex  hull  is  an  interval  on  the  line.  If  V  is  tire  intersection  of  two  UPPER  (or 
two  LOWER)  sides  that  determine  the  lines  y  =  m^x  +  b-j  and  y  =  rr^x  +  l>2  then  the 
transform  of  V  is  the  interval  between  (one-dimensional)  points  (m-j)  and  (m2).  This 
is  because  all  UPPER  (or  LOWER)  lines  of  support  at  V  must  have  a  slope  between 
(m-|)  and  (in o).  If  V  is  a  leftmost  or  rightmost  point,  then  the  set  of  slopes  of  the 
UPPER  (or  I.OWER)  lines  of  support  at  V  is  an  infinite  interval  on  the  line.  For 
example,  in  Figure  5-4  the  UPPER  transform  of  vertex  E  is  the  interval  (-co,e]  and 
the  LOWER  transform  of  vertex  E  is  the  interval  [d/o),  where  "e"  is  the  slope  of 
side  e  and  "d"  is  the  slope  of  side  d. 

The  transform  gives  us  all  the  information  we  need  to  generate  the  0(N)  pairs  of 
antipodal  vertices.  For  example,  in  Figure  5-4,  if  L 1  is  the  line  determined  by  side  a, 
then  the  parallel  line  of  support  L2  passes  through  the  vertex  F.  Equivalently,  the 
transform  maps  side  a  to  point  a  and  vertex  F  to  interval  F  such  that  point  a  lies 
inside  the  interval  F.  Since  side  a  is  hounded  by  vertices  A  and  B,  we  have  the 
antipodal  pairs  of  vortices  (A,F)  and  (fi,F).  We  can  generate  all  0(N)  pairs  of 
antipodal  vertices  by  finding  which  intervals  contain  the  N  points  a.  b,  c,  etc. 

Theorem  ?fi:  Given  a  convex  polygon  of  N  sides,  we  can  generate  the 
0(*N)  pairs  of  antipodal  vertices  in  0(M)  time. 

Proof:  When  we  generate  the  UPPER  and  LOWER  sets  of  points  on  the 
line  (Equation  14),  they  will  ho  in  sorted  order  because  the  slopes  of  the 
sides  of  the  convex  hull  are  already  sorted.  We  can  thus  easily  scan  the 
two  sets  to  determine  which  interval  each  point  lies  in  and  ultimately 
generate  the  0(N)  antipodal  pairs  of  vertices  in  0(N)  time.  □ 

In  summary,  v/e  hove 

Theorem  27:  Wc  can  compute  the  diameter  of  N  planar  points  in 
0(N  log  N)  worst-case  time. 

Proof:  We  first  construct  the  convex  hull  in  0(N  log  N)  ([43])  time  and 
then  compare  the  0(NI)  antipodal  pairs  of  vertices  in  0(N)  time.  To 
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generate  the  antipodal  pairs  we  may  use  either  the  scan  around  the 
convex  hull  of  Shames  [01]  or  first  transform  the  sides  of  tire  convex  hull 
to  points  on  a  line  and  then  perform  an  equivalent  scan  of  those  points. 

□ 

5.2.2.  Diameter  in  Three  Dimensions 

Many  features  of  our  algorithm  for  the  diameter  of  a  two-dimensional  set  of  points 
extend  to  three  dimensions.  We  first  construct  the  convex  hull  of  the  N  points  in 
0(N  log  N)  time  (by  the  algorithm  of  Preparata  and  Hong  [S3])  to  enable  us  to  solve 
the  diameter  of  the  set  of  points  as  the  diameter  of  a  convex  hull.  To  find  the 
diameter  of  the  convex  hull  we  then  generate  tire  set  of  antipodal  pairs  of  vertices, 
one  pair  of  which  determines  the  diameter.  In  this  section  we  show  how  to  generate 
the  K  pairs  of  antipodal  vertices  in  0((N  +  K)  log  N)  time  and  that  we  can  thereby 
compute  the  diameter  of  N  points  in  three-space  in  0((N  +  K)  log  N)  time. 

In  the  plane  we  used  the  concept  of  line  of  support  to  generate  the  O(N)  pairs  of 
antipodal  vertices.  In  three  dimensions  the  corresponding  concept  is  plane  of 
support.  For  each  face  of  the  convex  hull,  say  face  PQR,  the  plane  of  suppor 
passes  through  a  vertex  W  of  the  convex  hull.  This  generates  the  antipodal  pairs 
of  vertices  (P,W),  (Q.W),  and  (R,W).  Although  the  convex  hull  has  only  O(N)  faces,  in 
the  worst  case  there  may  still  he  pairs  of  antipodal  vertices.18  When  the 

number  of  pairs  K  is  less  than  Q(N^),  though,  it  is  not  obvious  how  to  generate  them 
in  less  than  OfN'’)  time. 

Our  first  step  toward  generating  the  antipodal  vertices  is  to  divide  the  faces  of 
the  convex  hull  into  the  two  sets  UPPER  and  LOWER.  (The  plane  that  contains  an 
UPPER  face  lies  above  the  convex  hull  and  the  plane  that  contains  a  LOWER  face 
lies  below  the  convex  hull.)  This  division  has  the  property  that  if  a  plane  L  contains 
a  face  of  the  UPPER  set,  then  the  plane  of  support  M  that  is  parallel  to  L  passes 
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through  a  vortex  of  the  i  OWE  It  set.  Similarly,  n  face  of  the  LOWER  set  determines  a 
parallel  plane  of  support  that  passes  through  a  vertex  of  the  UPPER  set.  (We 
consider  the  vertices  on  the  hnuiulary  of  the  UPPER  and  LOWER  sets  to  belong  to 
both  sets.)  We  next  show  how  wc  can  transform  the  UPPER  and  LOWER  sets  to  find 
the  antipodal  vertices  quickly. 

5.2. 2.1.  Transform  in  Three  Dimensions 

To  generate  the  0(N)  pairs  of  antipodal  vertices  we  use  an  extension  of  the 

geometric  transform  that  we  used  for  the  two-dimensional  problem.  In  two 

dimensions  our  choice  of  transform  was  motivated  by  the  fact  that  the  search  for 

lines  of  support  involves  only  comparisons  of  the  slopes  of  the  sides  of  the  convex 

hull.  Similarly,  in  three  dimensions,  the  search  for  planes  of  support  involves 

comparisons  of  the  slopes  of  the  faces  of  the  convex  hull.  Our  transform  for  the 

faces,  edges,  and  vertices  of  the  convex  hull  ail  follow  the  same  schema: 

The  UPPER  transform  of  a  component  of  the  convex  hull  is  the  set  of 
slope-pairs  of  all  UPPER  planes  of  support  that  contain  it.  Similarly,  the 
LOWER  transform  is  the  set  of  slope-pairs  of  all  LOWER  planes  of  support 
that  contain  it. 

Since  only  one  plane  of  support  contains  a  face  of  the  convex  hull,  a  face  maps  to 
one  point.  An  edge  is  contained  by  a  set  of  planes  of  support  with  one  degree  of 
freedom  so  an  edge  maps  to  an  interval  of  a  line.  Finally,  a  vertex  of  the  convex 
hull  is  contained  by  a  set  of  planes  of  support  with  two  degrees  of  freedom  so  a 
vertex  maps  to  a  planar  region.  We  now  present  algebraic  descriptions  of  the 
transforms  of  faces,  edges,  and  vertices  of  the  convex  hull. 

Transform  of  a  Face 

Let  z  -  ax  +  by  +  o  ho  the.  plane  determined  by  a  face  F  of  the  convex  hull.  The 
transform  of  face  F  is  simply  the  pair  of  slopes  (a,b),  a  point  in  the  ab  plane.  We 
may  alternately  view  this  transform  ns  a  combination  of  the  point  /'flat  duality  and 
orthographic  projection: 

2  =  ox  +  by  +  c  -*  (o,b,c)  -* 


(a.b). 


(15) 
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We  next  describe  the  transform  of  an  edge  and  a  vertex.  Since  the  UPPER  and 
LOWER  transforms  arc  so  similar,  wc  present  only  the  UPPER  transform  'for  each 
case. 

UPPER  Transform  of  an  Edge 

Suppose  that  two  faces  L  and  M  of  the  convex  hull  have  an  UPPER  edge  Ey^  in 
common.  Let  the  two  planes  L  and  M  determined  by  these  faces  be  written 


z  =  aLx  +  bLy  +  cL,  and 

z  s  aux  +  bv,y  ♦  <V 

We  can  write  the  line  where  planes  L  and  M  meet  in  parametric  form  as 


(16) 


(Px.Py,Ps)  ♦  u 


I  J  K 

aL  bL  -1 
®M  bM 


u  f  Reals 


(17) 


where  P  a  (Px,py,pz^  is  son,e  point  in  botb  P,an®s  ar*d  L  X  and  K  are  unit  vectors 
parallel  to  the  x,  y,  and  z  axes,  respectively.  For  example,  if  the  line  intersects  the 
xy-plane  we  can  choose  P  as 


P  *  (Px.Py.Pc)  * 
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A  plane  z  *  ax  ♦  by  ♦  c  contains  the  line  of  Equation  (17)  iff 

a  b  -1 
aL  bL  -1 

°M  bM  “1 

Since  the  transform  of  Eyyj  is  limited  to  UPPER  planes  of  support  that  contain  Elm. 
we  must  restrict  the  solutions  of  Equation  (13)  to  an  interval  defined  by  the  two 
points  (ay  by.  c.L)  and  (nM,  hM,  cM).  The  transform  of  edge  Ej_m  <*  the  projection  of 


-  0  and  P. » 


aP„  ♦  bPy  ♦  c. 


(18) 
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art  interval  of  this  line  to  the  xy  (ah)  plane.  We  summarize  our  results  in  the 
following  theorem: 

Theorem  PC:  Let  the  planes  L  and  M  of  Equation  (1G)  be  determined  by 
two  faces  of  a  convex  polyhedron  that  meet  at  an  UPPER  edge  Eyyj.  The 
UPPER  transform  of  edge  E^  is  the  set  of  points  (a,b)  on  the  line 

(bM  -  bL)a  -  (aM  -  aL)b  =  aLbM  -  aMbL  (19) 

that  lie  in  an  interval  determined  by  the  points  (a|_,  b[_)  and  (n^,  b^).  If  L 
and  M  arc  both  UPPER  planes  then  the  interval  lies  between  the  two 
points.  If  L  is  UPPER  but  M  is  LOWER  then  the  interval  is  the  ray  from 
point  (a|_,l>L)  that  docs  not  include  point  (a^.b.%1).  (Similarly  when  L  is 
LOWER  and  M  is  UPPER.) 

Proof:  The  line  of  Equation  (10)  is  from  Equation  (IS).  Since  the  set  of 
UPPER  planes  of  support  that  contain  edge  E^  is  connected,  the  UPPER 
transform  of  edge  is  an  interval  of  this  line.  Since  the  planes  L  and 
M  of  Equation  (10)  contain  edge  E^.  the  points  (oj_.bj_)  and  (a^.b^)  must 
be  on  the  UPPER  transform  of  edge  E^  iff  planes  L  and  M  are  UPPER 
faces  of  the  convex  hull.  Furthermore,  since  the  faces  L  and  M  are  the 
extreme  limits  that  a  plane  of  support  can  be  rotated  about  edge  Ey^, 
point  (aL,l»L)  must  be  an  endpoint  of  the  interval  if  plane  L  is  an  UPPER 
plane  of  support  and  (f|y;.l>M)  must  be  an  endpoint  of  the  interval  if  plane 
M  is  an  UPPER  plane  of  support.  From  these  conditions  it  follows  that  if 
both  L  and  M  are  UPPER  planes  of  support,  then  the  interval  must  be 
between  the  two  points.  If  L  is  an  UPPER  plane  of  support  but  not  M, 
then  the  interval  is  a  ray  starting  from  point  (oL.bL)  that  does  not  contain 
point  (Of^.bf^).  Similury  if  L  is  LOWER  and  M  is  UPPCR.  □ 

UPPER  Transform  of  a  Vertex 

Suppose  that  II  faces  meet  at  a  vertex  V  =  (Vx.Vy,Vz)  on  the  UPPER  part  of  the 
convex  hull.  Number  these  faces  in  counterclockwise  order  so  that  the  edges  that 
meet  at  vertex  V  are  l'^,  ,  . . .,  Pj.j -j .  (If  vertex  V  lies  on  the  boundary  of  the 

UPPER  and  LOWER  parts  then  include  only  those  edges  in  the  UPPER  part.)  Let  the 
planes  determined  by  the  It  faces  bo. 

z  =  OjX  +  bjy  +  Cj,  i  =  1 . H. 

If  ray  rj  is  the  ray  that  originates  at  point  V  and  points  down  edge  Ey+i  then 

ri  s  <VX.  V  V  ♦  u  fej.  ft.  7i).  u  >  0 


(20) 
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where 

*i  bi  -1 

a,*i  bi*i  -1  •  *21) 

ai+2  b.+2 

(Note  th.it  the  subscripts  "i+1"  and  "i+2"  are  to  be  taken  modulo  H.)  The  first 
determinant  in  equation  (21)  determines  the  line  that  ray  r-,  lies  in  and  the  second 
determinant  determines  which  of  the  two  possible  directions  r;  should  point.  We  may 
now  characterize  the  set  of  UPPER  planes  of  support  of  UPPER  vertex  V. 

Theorem  20:  Let  V  =  (Vx.Vy,VJ!)  be  an  UPPER  vertex  of  a  convex 
polyhedron  at  which  the  UPPER  rays  rj  determined  by  V  satisfy  Equation 
(20).  A  plane  z  *  ax  ♦  by  ♦  c  is  an  UPPER  plane  of  support  at  vertex  V  iff 

V2  *  a Vx  +  bVy  ♦  c  (22) 

and 

7j  £  aaj  +  bjjj.  i  »  1 . H.  (23) 

Proof:  A  plane  z  *  ax  ♦  by  ♦  c  is  an  UPPER  plane  of  support  at  vertex  V 
iff  it  passes  through  V  and  remains  above  all  the  rays  rj.  Equation  (22) 
requires  the  plane  to  pass  through  V.  The  plane  remains  above  the  H  rays 
iff 

Vz  ♦  uyj  £  a  (V.x  +  u^)  +  b  (Vy  ♦  Ujjj)  ♦  c,  Vu  >  0.  (24) 

By  subtracting  Equation  (22)  and  dividing  by  u  we  obtain  the  inequality 
(23).  Conversely,  if  Equations  (22)  and  (23)  ore  satisfied,  we  easily 
derive  Equation  (24)  by  multiplying  (23)  by  u  £  0  and  adding  to  (22).  □ 

For  the  transform  of  UPPER  vertex  V  vve  arc  interested  only  in  the  slopes  of  the 
UPPER  planes  of  support  at  V.  By  Theorem  29  we  have 

Theorem  30:  Let  V  be  an  UPPER  vertex  of  a  convex  polyhedron  at 
which  the  UPPER  rays  r4  determined  by  V  satisfy  Equation  (20).  The 
UPPER  transform  of  V  is  a  convex  polygonal  region  of  the  ab  plane 
determined  by  the  inequalities 

Yj  £  a«j  ♦  bflj,  i  *  1 . H. 


I  J  K 

(ec,.0,.r .)  s  +  +  s  ai  bi  *1  *  sgn 

ai*i  bi+i 
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5. 2. 2. 2.  Algorithm  for  Generating  Antipodal  Vertices 

We  have  just  seen  how  to  transform  the  faces  of  a  convex  polyhedron  to  points 
in  the  plane,  edges  to  edges  in  the  plane,  and  vertices  to  convex  polygonal  regions. 
Combining  these  we  can  transform  both  the  UPPER  and  LOWER  parts  of  the  convex 
hull  to  outerplanar  straight-line  graphs  in  the  plane,  as  illustrated  in  Figure  5-5.  We 
now  present  the  property  of  these  graphs  that  allows  us  to  find  efficiently  the  K 
pairs  of  antipodal  vertices. 

Let'V  be  on  UPPER  vertex  of  a  convex  polyhedron.  The  UPPER  transform  of  V  is  a 
convex  polygonal  region  V'  of  the  ah  plane.  By  definition  of  V',  any  UPPER  plane  that 
passes  through  vertex  V  transforms  to  a  point  in  V'  iff  it  is  an  UPPER  plane  of 
support.  Any  LOWER  plane  of  support  that  maps  to  a  point  in  V*  is  therefore  a 
parallel  plane  of  support  for  some  UPPER  plane  of  support  that  passes  through  V.  If 
a  LOWER  vertex  W  maps  to  a  region  W*  such  that  V*  and  W*  overlap,  then  V  and  W 
share  parallel  planes  of  support  and  are  therefore  antipodal. 


Figure  5-5:  Search  for  planes  of  support  maps  to  locating  overlapping  regions. 

We  hove  reduced  the  problem  of  computing  the  diameter  of  N  points  in 
three-space  to  the  problem  of  finding  intersections  of  the  regions  of  two  outerplanar 
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straight-line  graphs.  We  thus  have  the  following  outline  for  a  three-dimensional 
diameter  algorithm: 

Outline  of  Thrcc-Diracnsionat  Diameter  Algorithm 

1 .  Construct  the  convex  hull  of  the  N  points. 

2.  Divide  the  convex  hull  into  the  two  parts  UPPER  and  LOWER. 

3.  Map  both  the  UPPER  and  LOWER  parts  to  outcrpl&nar  straight-line 
graphs  by  the  transform  described  above. 

4.  Find  all  of  the  pairs  of  regions  in  the  UPPER  and  LOWER  graphs  that 
overlap.  Update  the  maximum  distance  between  pairs  of  antipodal 
vertices. 

5.  Report  the  maximum  distance  measured  as  the  diameter. 

We  can  construct  the  convex  hull  (Step  1}  in  0(N  log  N)  time  [83]  and  in  O(N) 
time  partition  the  hull  into  the  sets  UPPER  and  LOWER  (Step  2)  and  transform  these 
sets  to  outcrplanar  straight-line  graphs  (Step  3).  Since  there  are  only  O(N)  pairs  of 
antipodal  vertices.  Step  5  requires  only  constant  time.  This  leaves  Step  4,  finding 
the  overlaps  between  the  UPPER  and  LOWER  outcrplanar  graphs.  We  now  describe 
how  to  conduct  this  search  in  0((N  +  K)  log  N)  time,  where  K  is  the  number  of  pairs 
of  antipodal  vertices. 

Assume  that  there  are  0(N)  faces  on  both  the  UPPER  and  LOWER  parts  of  the 
convex  hull.  The  suhproblcm  of  Step  4  is  thus  to  find  the  overlaps  among  two  sets 
of  O(N)  convex  regions  (with  a  total  of  0(N)  edges).  Two  regions  overlap  iff  (1)  a 
vertex  of  one  region  is  contained  in  the  other  region  or  (2)  an  edge  of  one  region 
intersects  an  edge  of  another  region.  There  exist  algorithms  that  solve  these  two 
cases  separately;  we  only  need  to  combine  them. 

We  can  solve  the  first  case,  determing  inclusion  of  the  vertices  in  planar  regions, 
by  a  number  of  algorithms.  Lee  and  Preparata  [GG]  describe  how  we  can  locate  a 
point  in  the  correct  region  in  O(iog^N)  time,  with  only  0(N  log  N)  preprocessing  time. 
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We  can  thus  locate  the  0(N)  points  in  0(N  locj^N)  time.  Lipton  and  Tarjan  [70]  have 
improved  the  query  time  to  0(log  N),  yielding  an  0(N  log  N)  time  algorithm  for  the 
O(N)  points  while  using  only  0(N)  storage,  but  the  "constant  factor"  of  their 
algorithm  is  very  large.  Preparata  [SO]  has  produced  a  practical  algorithm  that 
costs  only  6  f  log  N  1  comparisons  for  each  query  but  requires  0(N  log  N)  storage 
and  preprocessing  time.  Since  we  want  to  locate  a  set  of  points  together  rather 
than  just  one  at  a  time,  however,  we  can  use  Preparata's  0(K  log  K)  + 
O(N)  +  0(K  log  N)  time  algorithm  [82]  or  Lee  and  Yang's  0((N  +  K)  log(N  +  K))  time 
algorithm  [68]  for  locating  a  set  of  K  points  in  a  straight-line  planar  graph  of  N 
vertices,  giving  us  an  0(N  log  N)  time  algorithm  for  locating  O(N)  points. 

The  second  case,  finding  all  K  intersections  of  the  edges  of  the  two  graphs,  can 
be  solved  in  0((N  +  K)  log  N)  time  and  0(N)  storage  by  Brown's  modification  [24]  of 
Bentley  and  Ottmonn's  [12]  intersection  algorithm  (which  is  itself  a  modification  of 
Hoey's  algorithm  for  determining  if  any  of  N  line  segments  intersect  [94,  91]). 

The  algorithm  below  finds  all  of  the  edge  intersections  and  all  of  the  vertex 
inclusions  in  0((N  +  K)  log  N)  time  and  0(N)  storage.  It  finds  the  edge  intersections 
by  Brown's  algorithm  [24]  and  the  regions  containing  the  vertices  by  maintaining  two 
order  relations  (Ry  and  R|_)  for  the  two  outerplonar  graphs.  Since  the  extra  time 
required  to  maintain  Ry  and  Ry  is  only  0(N  log  N)  and  the  storage  is  only  O(N),  the 
total  time  and  storage  bounds  are  the  same  as  for  the  edge  intersection  algorithm. 
This  algorithm  uses  several  simple  data  structures  and  functions: 
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-  Gy  and  Gy  -  two  outerplanar  straight-line  graphs  representing  the 
transforms  of  the  UPPER  and  LOWER  parts  of  the  convex  hull  of  the 
three-dimensional  points, 

>  Nextlntfi],  »  next  detected  intersection  point  for  segment  i  (that  is  to 
the  right  of  the  current  x  coordinate  of  the  !eft-to*right  scan), 

-  Q  -  a  queue  of  (<intersection-point  or  endpoint>,  <segment>)  pairs, 
sorted  by  the  x  coordinates  of  the  intersection  points  (or  endpoints), 

-  Ry  [Ry]  -  an  order  relation  of  edges  from  Gy  [Gy]  (evaluated  at  the 
current  x  coordinate  of  the  left-to-right  scan), 

-  R  -  an  order  relation  of  line  segments  for  both  Gy  and  Gy  combined 
(evaluated  at  the  current  x  coordinate  of  the  left-to-right  scan), 

.  / 

-  Regions[i]  «  ordered  pair  (Rogions[i].above,  Regions[i].be!ow)  of  the 
regions  above  and  below  edge  i  in  the  graph  (Gy  or  Gy)  of  edge  i, 

-  Point[i]  *  (Xj,  yjt  Zj)  of  the  three-dimensional  point  corresponding  to 
region  i  of  Gy  or  Gy, 

-  Vertices[i]  *  (three)  vertices  of  the  convex  hull  that  determine  the 
face  that  maps  to  vertex  i  of  Gy  or  Gy. 

-  lnsert(P,A,G)  -  inserts  (P,A)  into  the  queue  Q,  where  P  is  the  left  or 
right  endpoint  of  segment  A  or  the  intersection  of  segment  A  and 
another  segment. 

-  Delete(P,A,Q)  -  deletes  (P,A)  from  the  queue  Q.  Exception:  If  P.x  *  co 
the  request  is  ignored. 

-  Insert(S.R)  [Oelete(S,R)]  -  inserts  [deletes]  segment  S  in  [from]  order 
relation  R  (where  the  order  is  evaluated  at  the  x  coordinate  of  the  left 
endpoint  of  S). 

-  Above(P.R)  [Below(P.R)]  -  returns  the  segment  above  [below]  point  P 
in  order  relation  R  (where  the  order  Is  evaluated  at  P.x), 
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-  lnsert(S.Ry,R|_)  [Delete(S,Ry.Ry)]  -  inserts  [deletes]  segment  S  in 
[from]  the  order  relation  (Ry  or  R|_)  to  which  S  belongs  (where  the 
order  is  evaluated  at  the  x  coordinate  of  the  left  endpoint  of  S), 

-  Above(P,Ry,RL)  [Below(P,r.y,R|_)]  *  returns  the  segment  above  [below] 
point  P  in  the  order  relation  (Ry  or  R[_)  to  which  P  does  not  belong 
(where  the  order  is  evaluated  at  P.x). 

-  Pairs(L,M,D)  -  computes  the  distances  between  all  pairs  of  points  (i,j) 
such  that  i  C  L  and  j  {  M.  If  any  of  these  distances  are  greater  than  D, 
then  D  is  set  to  the  new  maximum. 


Algorithm  for  Search  Step  of  Diameter  Algorithm 


proc  Inter(A.B) 

!  Implement  the  modified  insertion  rule  of  [23]  for  segments  A  and  B.; 
P  ♦-  Intersec lio»(A,B); 

L  *■  Nextlnt[A];  M  <-  Nextlnt[B]; 
if  P.x  <  L.x  then 

Delete(L,A,0);  lnsert(P,A,Q);  NextlntfA]  P; 
if  P.x  <  M.x  then 

Delete(M,B,Q);  lnsert(P,B,Q);  Nextlnt[B]  <-  P; 


!  Initialization; 

Q  «-  {  all  pairs  (P,i)  where  P  is  a  left  or  right  endpoint  of  segment  i, 
sorted  by  the  x  coordinates  P.x  } 

R  <-  Ry  «-  Ry  <-  <t>;  !  order  relations  for  the  edges  of  Gy  and  G(_; 

Nextlnt[i].x  <•  co  for  all  segments  i; 

Diam  <-  0; 
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while  Q  j*  <t>  do 

(P,S) «-  (next  (point, segment)  pair  on  Q);  !  X  coord  of  scan  becomes  P.x; 

T  «-  (other  segment  intersecting  at  P  If  the  next  pair  on  Q  is  (P,T)); 

If  P  is  the  left  endpoint  of  segment  S  then 
Insert(S.R); 

A  «-  Above(P,R);  B  «-  BalovHP.R); 

If  A  intersects  S  then  Inter(A.S); 

If  B  intersects  S  then  Inter(B.S); 

Pair3(Point[Regions[Above(P1Ru,P.L)].below],  Vertices[P],  Diam);  !  Three  pairs; 

Else  if  P  is  the  right  endpoint  of  segment  S  then 
A  ♦-  Abovc(P.R);  B  ♦-  Bciow(P,R); 

Delete(S.R); 

If  A  intersects  B  then  Inter(A.B); 

Pairs(Point[RRgions[Above(P,Ry,P.y)].below],  Vcrtices[P],  Diam);  !  Three  pairs; 

Else  !  P  is  an  intersection  of  segments  S  and  T.; 

Report(P); 

Nextlnt[S].x  Nextlnt[T].x  «-  co; 

Reverse(S,T,R);  I  Let  S  become  the  top  segment.; 

A  «-  Above(P.R);  B  4-  Below(P.R); 

If  A  intersects  S  then  Inter(A.S); 

If  B  intersects  T  then  Inter(B.T); 

Pairs(Regions[S],  Regions[T],  Diam);  !  Compare  four  pairs.; 

We  have  just  seen  how  to  find  all  of  the  0(K)  overlaps  of  regions  of  two 

outerplanar  straight-line  graphs  of  size  0(N)  in  0((N  +  K)  log  N)  time  and  O(N) 

storage.  But,  as  shown  above,  this  implies 

Theorem  3 1 :  The  diameter  of  N  points  in  three-space  can  be  computed 
in  0((N  ♦  K.)  log  N)  time,  where  K  is  the  number  of  pairs  of  antipodal 
vertices  on  the  convex  hull  of  the  N  points. 


5.2.3.  Refinements,  Extensions,  Related  and  Unsolved  Problems 

We  have  concentrated  only  on  worst-ease  two-  and  three-dimensional  algorithms 
for  computing  the  Euclidean  diameter  exactly.  In  this  section  we  briefly  describe 
results  for  fast  expected-time,  approximation,  and  higher-dimensional  algorithms, 
open  problems  and  an  application  to  Chebyshev  regression. 
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1 .  One  of  the  major  open  problems  for  the  Euclidean  diameter  is  proving  a 
nontrivial  lower  bound.  In  Appendix  II  we  show  that  a  diameter 
algorithm  can  solve  an  empty-intersection  problem  for  which  an 
Q(N  log  N)  time  lower  bound  has  been  proven  for  a  weak  model  of 
computation,  but  there  is  still  no  Q(N  log  N)  time  lower  bound  for  a 
model  of  computation  strong  enough  to  construct  a  convex  hull  in 
0(N  log  N)  time.  Shamos  [92]  conjectures  that  the  diameter  problem 
has  a  worst-case  lower  bound  of  o(N  log  N)  time  for  any  metric  whose 
circle  has  a  continuously-turning  tangent  (such  as  the  Euclidean 
metric),  but  that  if  the  circle  has  only  a  discretely-turning  tangent 
(such  as  the  L-j  or  Laj  metrics)  then  we  can  compute  the  diameter  in 
O(N)  time. 

2.  In  our  two-dimensional  diameter  algorithm  the  most  expensive  step  is 
the  construction  of  the  convex  hull  of  the  N  points  in  0(N  log  N) 
worst-case  time.  A  fast  eixpactoxl-tima  convex  hull  algorithm  leads  to 
a  fast  diameter  algorithm  [91].  For  example,  if  the  expected  number 
of  vertices  on  the  convex  hull  is  only  0(NP)  for  some  p  <  1,  then  we 
may  construct  the  convex  hull  in  0(N)  expected-time  [16]  and 
therefore  compute  the  diameter  of  N  points  in  0(N)  expected-time. 

Our  three-dimensional  algorithm,  however,  takes  not  only  0(N  log  N) 
worst-case  time  to  construct  the  convex  hull,  but  also  0((N  +  K)  log  N) 
time  to  generate  the  K  pairs  of  antipodal  vertices.  If  the  expected 
number  of  vertices  on  the  convex  hull  is  only  0(NP)  for  some  p  <  1  and 
the  expected  number  of  pairs  of  antipodal  vertices  is  K,  then  we  may 
compute  the  diameter  in  0(N  +  K  log  N)  expected-time.  An  obvious 
open  problem  is  to  prove  bounds  for  K,  as  a  function  of  N  for 
interesting  distributions  of  points.  (Conjecture:  Let  the  straight-line 
planar  graphs  that  the  UPPER  and  LOWER  parts  of  the  convex  hull  map 
to  be  Gy  end  G|_  and  let  the  regions  of  Gy  be  Uj  and  the  regions  of  Gy 
be  Lj.  Fot  any  hounded  region  m  of  Gy  or  Gy  let  R(m,0)  be  the  aspect 
ratio  of  m  —  the  ratio  of  height  to  width  —  when  m  is  rotated  an  angle 
0.  For  any  unbounded  region  m  with  sides  along  rays  r  and  s,  let 
R(m,0)  be  the  aspect  ratio  for  any  (rotated)  isosceles  triangle  with  its 
two  equal  sides  on  lines  r  and  s.  The  number  of  overlapping  regions  K 
is  O(S^^N),  where  S  =  max(i,j,0)  R(Uj,Q)  /  R(Lj,0).  This  type  of  bound 
arises  in  the  maximum  overlap  of  N  rectangular  regions  [22].) 

In  D  dimensions  we  can  construct  the  convex  hull  of  N  points  in  0(N) 
expected-lime  if  the  D  coordinates  of  the  points  are  independently 
distributed  [10]  [30].  In  this  case  the  expected  number  of  vertices 
on  the  convex  hull  is  only  0(locp_1N)  [10]  and  we  can  then  compute 
the  diameter  in  only  OOog*^'1  ^N)  more  expccted-time  by  the 
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brute-force  method  of  comparing  all  pairs  of  vertices. 

3.  The  best  known  worst-case  three-dimensional  Euclidean  diameter 
algorithm  is  that  of  Yao  [102],  which  runs  in  OfN1’®)  time.  Yao  has 
also  produced  a  D-dimensional  Euclidean  diameter  algorithm  that  runs  in 
0(n2-«(D))  time  where  a(D)  s  2"^+^>.  (This  algorithm  does  not 
involve  construction  of  a  convex  hull  because  a  convex  hull  in  four  or 
more  dimensions  may  have  0{N^)  edges  ([49],  p.193).) 

4.  Yuvnl's  0(N  log  N)  time  two-dimensional  Chebyshev  regression 
algorithm  [01]  relies  on  a  scan  of  a  convex  hull  with  lines  of  support. 
We  may  apply  our  0((N  +  K)  log  N)  time  algorithm  for  locating  planes  of 
support  of  a  convex  polyhedron  to  obtain  an  0((N  +  K)  log  N)  time 
three-dimensional  Chebyshev  regression  algorithm. 

5.  We  may  approximate  the  Euclidean  diameter  of  a  planar  set  of  N  points 
within  a  factor  of  1  ♦  (  in  worst-case  time  0(M  +  I/O  [7]  (Section 
6.1).  (Shamos  and  Yuval  [95]  have  previously  described  a  technique 
that  leads  to  an  0(M/((1/2))  time  approximation  algorithm.)  Bentley, 
Faust  and  Preparata  [7]  describe  a  D-dimensional  (-approximate 
algorithm  that  runs  in  0(N  +  (1/6)2(D-1))  tjme  Open  Problem: 
Construct  a  faster  0-dimensional  (-approximate  diameter  algorithm  or 
prove  a  nontrivial  lower  bound  for  the  D-dimensional  (-approximate 
diameter  problem. 


€.3.  Summary 

For  both  of  the  problems  treated  in  this  chapter  (linear  programming  and  the 
Euclidean  diameter  of  a  set  of  points)  we  have  found  that  an  apparently 
D-dimensional  problem  can  be  expressed  as  a  D-1  “dimensional  problem.  This  is 
because  each  problem  can  be  expressed  as  a  search  for  a  flat  (or  a  pair  of  flats)  In 
which  only  the  0-1  slopes  (but  not  the  intercept)  are  important.  The  transform  is 
chosen  not  only  to  reduce  the  dimensionality  of  the  problem  but  also  to  represent  it 
in  a  form  for  which  there  are  already  fast  algorithms.  The  point  /  flat  duality 
transforms  the  search  of  flats  into  a  search  of  points  and  orthographic  projection 
removes  the  unneeded  coordinate. 
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6.  Miscellaneous  Problems  and  Techniques 

In  this  chapter  wn  present  two  problems  and  solutions  that  do  not  fit  into  any  of 
the  categories  of  the  previous  chapters.  The  techniques  applied  to  the  first 
problem  (finding  the  approximate  diameter  of  a  set  of  planar  points)  provide  a 
contrast  to  the  techniques  for  the  exact  diameter  in  the  previous  chapter  because 
the  emphasis  is  on  quickly  producing  an  approximate  convex  hull  rather  than  on  a 
fast  search  of  an  exact  convex  hull.  Our  solution  to  the  second  problem  (fitting 
points  on  a  hemisphere)  introduces  a  new  transform  called  gnomonic  projection  that 
has  the  property  of  mapping  great  circles  on  a  sphere  to  straight  lines  in  the 
Euclidean  plane. 

6.1.  Approximate  Diameter  of  Points  in  Two  Dimensions 

We  presented  an  0(M  log  N)  time  algorithm  for  computing  the  exact  Euclidean 
diameter  of  N  planar  points  in  Section  5.2.  In  this  .section  we  describe  two  O(N)  time 
algorithms  for  approximating  the  Euclidean  diameter  of  N  planar  points  within  a 
relative  factor  of  (1+#).19  To  achieve  the  faster  time  we  must  use  radically 
different  techniques.  Whereas  the  exact  diameter  algorithm  used  the  point  /  flat 
duality  and  orthographic  projection  to  obtain  a  problem  of  locating  points  within  a 
tesselation,  the  first  approximation  algorithm  uses  rotation  to  define  a  metric  whose 
unit  circle  is  a  regular  polygon  and  the  second  approximation  algorithm  extends  the 
ideas  of  the  first  by  defining  a  transform  based  on  a  "pie-slice"  diagram  and  use  of 
the  floor  function. 

6.1.1.  First  Approximate  Diameter  Algorithm 

Shamos  and  Yuvnl  [95]  described  how  to  approximate  the  Euclidean  metric  with  a 
metric  whose  unit  circle  is  a  regular  polygon  —  the  distance  between  two  points  P 
and  Q  is  the  width  of  the  smallest  regular  K-gon  (of  a  particular  orientation)  that 
contains  both  P  and  Q.  (Figure  G-1  illustrates  the  case  for  K  *  3  sides.)  This 


19 


'The  source  of  this  problem  a  content  held  by  *i.l.  Shamos  at  Carnegie- Mellon  University  Airing  Fall  1978% 
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approximation  has  tlio  properly  Hunt  as  the  number  of  sides  of  the  regular  polygon 
increases  the  approximation  improves.  Our  problem  is  to  find  a  function  K(0  such 
that  we  can  approximate  the  diameter  within  a  factor  of  1  +  i  by  using  a  regular 
polygon  metric  of  K(0  sides. 


Figure  6-1:  Points  A  through  H  determine  the  "octagon"  metric  diameter 


Theorem  32:  Let  DIAM^  be  the  diameter  of  a  planar  set  of  points  S  in  a 
regular  K-gon  metric  and  let  01AM  be  the  Euclidean  diameter  of  S.  If 

K  >  T  n  /  sec’ll  +  0  1 

then 


DIAMk  <  DIAM  <  (1  ♦  O  DIAMk. 


Proof:  The  worst  cases  are  pictured  in  Figures  6-2  (a)  and  (b).  In 
Figure  6-2  (a)  the  Euclidean  diameter  is  determined  by  points  P  and  Q  but 
the  K-gon  diameter  is  determined  by  points  T  and  U  and  points  V  and  W. 
Letting  D(A,R)  denote  the  Euclidean  distance  between  two  points  A  and  B, 
we  define 


r  =  D(T,U)  =  D(V,W)  and  R  =  D(P,Q). 
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Since  Uie  metric  is  based  on  ij  regular  Is-gon,  the  angle  od  between  lines 
Lj  and  Ly  and  between  lines  Ly  and  Lw  is 

od  *  2lt  /  K. 

The  minimum  possible  K-gon  diameter  r  is  achieved  when  points  T  and  U 
are  placed  so  that  line  TU  is  perpendicular  to  lines  Lj  and  Ly  and  points  \f 
and-W  arc  placed  so  that  line  VW  is  perpendicular  to  lines  Ly  and  Lyy. 
From  these  conditions  it  follows  that 

r  =  R  cos(cd/2). 

Since  (  ~  (R  -  r)  /  r,  we  have 

K  >  T  n  /  sec^O  ♦  O  1. 

In  Figure  6-2  (b)  the  Euclidean  diameter  is  determined  by  points  P  and 
Q  but  the  K-gon  diameter  is  determined  by  points  P  and  0  and  points  P 
and  C.  Letting  r  s  0(P,B)  =  0(P,C)  and  R  =  D(P,Q)  we  have 

«  .  in  .  .  (./r)...W/2)  -  1. 

r  r 

which  is  maximized  at  s/r  *  sin(od/2),  so  that 

«  S  cos(od/2)  +  sin(ed/2)  tan(od/2)  -  1  *  sec(od/2)  -  1 . 

Since  cd  ■  21T  /  K,  we  have 

K  >  fit  /sec’Ul  +Ol. 

□ 

Theorem  33^  If  K  s  f  n  /  sec'^O  ♦  O  1,  then  as  (  ■*  0, 
K  -*  Tt  /  (201/2. 

Proof:  Since  cos(x)  *  1  -  x2/2  +  . . .,  the  result  follows.  □ 

First  Approximate  Diameter  Algorithm 


1 .  Let  the  number  of  sides  of  the  regular  K-gon  metric  be 

K  *  T  TT  /  sec'll  ♦  20  1. 
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(A)  (B) 


Figure  6-2:  Worst  cases  for  K-gon  diameter. 

2.  For  each  of  the  directions  (angles)  0.  2n/K,  4rr/K,  6rr/K,  .  .  . 
2(K-1)n/K  find  the  most  extreme  of  the  N  points.  (This  is  equivalent 
to  repeatedly  rotating  the  N  points  and  finding  the  point  farthest  to 
the  right.) 

3.  Find  the  (exact)  diameter  of  the  (at  most)  K  points  determined  in  Step 
2  and  multiply  it  by  (1  +  ().  This  can  be  done  in  0(K)  time  by  a  Graham 
[48]  scan  to  construct  the  convex  hull  and  a  Shamos  [89]  scan  to 
compute  the  diameter. 

Theorem  34:  The  Euclidean  diameter  of  a  set  of  N  planar  points  can  be 
approximated  within  a  factor  of  (1  +0  in  0(N/(C-)1/::)  time. 

Proof:  By  Theorem  32  the  value  chosen  for  K  in  the  above  algorithm  is 
sufficient  for  the  (-approximation  of  the  diameter.  Theorem  33 
establishes  that  K  =  0(N/((),/3)*  Since  the  most  expensive  step  of  the 
algorithm  is  Step  2,  v/e  can  compute  the  (-approximate  diameter  in  0(KN) 
=  0(N/((),/2)  time.  □ 
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6.1.2.  Second  Approximate  Diameter  Algorithm 

Although  the  first  (-approximate  diameter  algorithm  takes  only  time  linear  in  N,  it 
is  also  linear  in  1/<l/?.  The  second  (-approximate  algorithm  reduces  the  time  from 

0(N/(()l/2)  to  0(N  +  1/()  by  using  a  transform  based  on  a  "pie-slice"  diagram  (Figure 
6-3)  and  the  floor  function.  (Bentley,  Faust,  and  Preparata  [7]  describe  a  different 
(-approximate  planar  diameter  algorithm  that  also  runs  in  time  0(N  +  I/O.)  Bentley, 
Weide,  and  Yao  [18]  have  used  a  simple  "pie-sHce"  diagram  for  their  Voronoi 
diagram  algorithm  and  Weide  [99]  has  used  the  floor  function  to  speed  up  some  of 
the  algorithms  in  his  thesis. 


Figure  6-3:  A  "pie-slice"  diagram. 

The  "pie-slice"  diagram  enables  us  to  find  efficiently  a  small  subset  S'  of  the  N 
points  of  S  such  that  the  diameter  of  S'  approximates  the  diameter  of  S.  The  center 
of  the  diagram,  where  the  K  "slices"  meet,  can  be  any  point  on  or  within  the  convex 
hull  of  S.  (The  slices  do  not  each  cover  an  angle  of  2TT/K  radians,  however, 
becatise  the  computation  of  which  slice  a  point  lies  in  would  then  require  inverse 
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trigonometric  functions.  The  slices  are  instead  chosen  to  divide  uniformly  the  slopes 
in  each  of  the  octants,  leading  to  simpler  compulations.)  From  each  of  the  slices  we 
choose  one  point  to  insert  in  the  subset  S'.  As  illustrated  in  Figure  G-3,  the  point 
that  we  choose  is  either  the  farthest  left,  farthest  right,  highest,  or  the  lowest, 
depending  upon  the  orientation  of  the  slice.  Wc  now  determine  how  large  K  must  be 
to  ensure  that  the  diameter  of  S'  is  within  a  factor  of  (1  +  0  of  the  diameter  of  S. 


(A)  (B) 


Figure  6-4:  The  two  cases  for  choice  of  diameter  with  the  "pie-slice"  algorithm. 


Theorem  39:  If  DIAM  is  the  diameter  of  a  set  S  of  N  planar  points,  then 
the  diameter  DIAM}-  of  the  points  chosen  from  a  pic-slice  diagram  of  K 
slices  satisfies 

DIAMk  <  DIAM  <  (1+ODIAMk 
if 


K  >  r  8  i/l  +fe/2  (1  +  /T+?)  /«!• 


Proof:  Tigure  0-4  illustrates  the  two  possible  worst  cases.  In  both 
cases  the  Euclidean  diameter  is  determined  by  points  P  and  Q  but  the 
approximation  algorithm  chooses  points  A  and  D  instead.  The  lines  AP  and 
BO,  however,  are  parallel  for  case  (A)  but  perpendicular  for  case  (B).  The 
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of  the  pie-slicc  dinoram. 

In  case  6-4  (A)  we  have 


P  =  (Px,  Py)  =  (Ax  -  <\(Ay-Cy),  Ay)  and 
Q  =  (Qx.Gy)  =  (Bx  -  \(By-Cy).  By). 


where 


X  =  8  /  K. 


Lettiny  DX  =  Ax  -  Bx  and  DY  =  Ay  -  By, 


x,2  D(P,Q)2  (DX  -  XDY)2  +  DY2 

(1  +  *)  =  — — -r  =  - 5 - - - 

D(A,B)2  DX2  +  DY2 


Equation  (25)  is  maximized  as  DY  /  DX  oo.  The  worst  case  for  Figure 
6-4  (A)  is  thus 

K  =  f  8  /  t1/:  1.  (26) 


For  case  6-4  (B)  we  have 


P  *  (Px,  Py)  *  (Ax.  Ay  -  X(Ax-Cx))<  and 

Q  *  (Qx.Qy)  =  (Bx  -  X(By-Cy).  By). 


(27) 


Although  the  coordinates  of  point  C  cancelled  out  for  case  (A),  case  (B)  is 
worst  when  the  points  A,  B,  and  C  form  an  equilateral  triangle  (Figure 
6-5).  We  now  have  only  to  determine  the  worst  orientation  for  this 
triangle.  If  C  is  the  origin  then  A  and  B  can  be  chosen  to  be  two  points  on 
the  unit  circle  centered  at  C ; 

A  -  (cos(8),  sin(0)),  and 
B  =  (ccs(0+n/3),  sin(G+n/3)). 


For  line  AP  to  remain  vertical  while  line  BO  remain.®  horizontal  (as  in  Figure 
6-4),  we  require  that 

Tt/4  >  0  *  -TT/12.  (29) 

Combining  Equations  (27)  and  (28)  we  obtain 
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(i  +<)2  =  =  d(p,Q)2 

D(A.B)2 

=  F  cos2(0)  +  G  sin2(0)  +  H  sin(0)cos(0), 


(30) 


where 

F  =  (4  +  6,/3\  *  7\2)  /  4, 

G  =  (4  +  2i/3\  +  X2)  /  4,  and 
H  =  (2\  +  yax2)  /  2. 

Equation  (30)  is  maximized  at 

tan(20)  =  3‘1/2, 
which  is  satisfied  in  the  range  of  Equation  (29)  at 

0  =  n/12.  (31) 

Plugging  Equation  (31)  into  Equation  (30)  we  have 

k  =  r  s  Vi  ♦  yf/2  (i  ♦  yr+T)  /  <  l  (32) 


Since  the  value  of  K  for  case  (A)  (Equation  (26))  is  only  linear  in  1/t1''2, 
the  worst  case  is  case  (B),  with  K  defined  by  Equation  (32).  □ 

With  a  little  calculation  we  have 

Theorem  36:  If  K  is  defined  by  Equation  (32),  then  as  €  -*  0, 


te'JS.  MMHK 
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Figure  6-5:  Worst  case  for  second  approximation  algorithm. 


Second  Approximate  Diameter  Algorithm 

1.  Pick  an  arbitrary  point  from  the  set  of  N  points  to  serve  as  the  center 
of  the  "pie-r.lice"  diagram.  (Actually,  any  point  within  or  on  the  convex 
hull  of  the  N  points  can  serve  as  the  center  of  the  diagram.)  For 
simplicity  of  exposition  we  will  assume  that  the  center  is  the  origin. 

2.  For  an  ^-approximation  to  the  Euclidean  diameter,  let  the  number  of 
regions  in  the  diagram  he 


K  =  f  4  ^/l  +^3/2  ( 1  +  1/TT2?)  /  (  1  . 
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3.  For  each  of  the  N  points  (xj.yj)  (a)  determine  which  bucket  j  the  point 
falls  in,  and  (b)  compare  Xj  (or  yj,  depending  on  the  bucket  j)  with  the 
most  extreme  value  yet  found  for  bucket  j.  We  can  easily  do  part  (b) 
in  constant  time,  but  to  avoid  a  binary  search  of  0(log  K)  steps  in  part 
(a)  we  must  use  the  floor  function.  For  the  first  octant  the  formula  is 

BUCKET  <-L  (K/S)  (y/x)  J  . 

For  the  other  octants  tire  formula  is  similar. 

4.  Find  the  exact  diameter  DIAM^'  of  the  K  extreme  points  (plus  the 
central  point  of  the  diagram)  determined  in  Step  3.  Using  a  Graham 
[48]  scan  to  construct  tire  convex  hull  and  then  a  Sliamos  [91]  scan 
to  find  the  hull's  diameter,  we  can  compute  DIAM^1  in  0(K)  time.  Return 

D!AMk  =  (1  +  <•)  DIAMk' 
as  the  estimated  diameter. 


Theorem  37:  An  approximation  DIAM^  of  the  Euclidean  diameter  of  N 
planar  points  that  satisfies 

DIAM  /  (1  +  O  <  D!AMk  <  (1+0  DIAM 
can  be  computed  in  0(N  +  1/0  time. 

Proof:  Theorem  35  defines  the  value  of  K  required  to  obtain  an 
(-approximate  diameter  from  a  "pie-slice"  diagram  and  Theorem  36  shows 
that  this  function  is  0(N  +  1/0.  □ 


6.2.  Fitting  Points  on  a  Hemisphere 

We  can  sometimes  solve  a  problem  involving  points  on  a  sphere  similarly  to  the 
corresponding  problem  involving  points  in  a  plane.  For  example,  the  spherical 
convex  hull  of  N  spherical  points  is  similar  to  a  planar  convex  hull  provided  that  the 
N  spherical  points  all  lie  within  a  hemisphere.  If,  however,  there  is  no  hemispherical 
cap  that  contains  all  N  points,  then  the  (interior  of  the)  spherical  convex  hull  is  the 
entire  sphere.  The  crucial  test  (for  the  spherical  convex  hull  problem)  is  to 
determine  if  there  exists  a  hemispherical  cap  that  contains  all  N  points.  This 
problem  is  a  special  case  of  the  problem  of  determining  the  densest  hemisphere 
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determined  by  N  points  on  a  K-dimensional  sphere.  Johnson  and  Preparata 
[56]  have  shown  tiic  densest  hemisphere  problem  to  be  NP-complete  when  N  and  K 
are  arbitrary.  (Tor  fixed  K,  however,  there  is  an  0(N*^  log  N)  time  algorithm.)  In 
this  section  wo  will  confine  ourselves  to  the  (simpler)  problem  of  determining  if  N 
spherical  points  can  all  be  fit  into  a  hemispherical  cap. 

There  are  three  0(N  log  N)  time  algorithms  for  determining  if  there  exists  a 
hemisphere  that  contains  a  given  set  of  N  spherical  points.  One  solution  is  to 
intersect  the  N  half-spaces  that  contain  the  sphere  and  whose  boundaries  are 
tangent  to  the  sphere  at  the  N  points.  If  the  resulting  polytope  is  unbounded,  then 
the  N  points  can  be  fit  into  a  hemispherical  cap.  A  second  solution  is  to  construct 
the  dual  of  the  spherical  Voronoi  diagram  in  0(N  log  N)  time  (Section  4.2)  and  test  if 
one  of  the  faces  of  the  diagram  can  contain  a  hemisphere.  But  the  most  interesting 
solution  uses  a  transform  called  a  gnomonic  projection  and  is  due  to  Yuval 
[105].  Yuvnl's  algorithm  uses  gnomonic  projection  to  transform  a  problem  that 
involves  great  circles  on  a  sphere  to  a  problem  that  involves  straight  lines  in  the 
Euclidean  plane.  To  start  simply,  we  will  first  solve  the  two-dimensional  case  and 
then  extend  our  result  to  three  and  four  dimensions. 

6.2.1.  The  Two-Dimensional  Caso 

The  two-dimensional  problem  is  to  determine  if  there  exists  a  semicircle  that 
contains  a  given  set  of  N  circular  points.  There  arc  two  ways  in  which  we  can  solve 
this  in  O(N)  time.  One  way  is  to  find  the  largest  gap  between  any  two  consecutive 
points  of  the  circle.  If  this  gap  is  greater  than  or  equal  to  a  semicircle,  then  all  N 
points  can  fit  in  the  other  semicircle.  Gonzalez  [47]  describes  how  to  find  the 
largest  gap  between  N  unsorted  points  on  a  line  in  0(N)  time  (by  using  the  floor 
function).  With  a  small  modification,  his  algorithm  will  find  the  largest  gap  on  a  circle 
and  thus  solve  the  two-dimensional  problem  in  0(N)  time. 

The  largest  gap  solution  does  not,  however,  extend  easily  to  the 
three-dimensional  (spherical)  problem,  whereas  the  gnomonic  projection  solution 
does.  Figure  G-6  illustrates  how  the  gnomonic  projection  algorithm  works.  Let  C  be 
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Figure  6-6:  Gnomonic  projection  of  points  on  a  circle. 


the  center  of  the  circle  ami  let  L  be  a  line  tangent  to  the  circle.  The  gnomonic 
projection  transforms  each  circular  point  P  to  a  point  P'  on  the  tangent  line  L  such 
that  points  P,  P',  and  C  are  collincar.  This  transform  is  not  one-to-one  because  for 
each  point  P'  on  line  L,  there  are  two  points  on  the  circle  that  map  to  P*.  We  can 
alleviate  this  ambiguity  by  giving  each  o'  the  N  projected  points  on  L  one  of  the  two 
labels  "red"  or  "green".  If  a  circular  point  P  is  on  the  top  half  of  the  circle,  then  it 
maps  to  a  point  P'  labelled  red.  If  point  P  is  on  the  bottom  half  of  the  circle,  then  it 
maps  to  a  point  P'  labelled  green.  We  now  use  gnomonic  projection  to  transform  our 
circular  problem  to  a  linear  problem  that  we  can  solve  in  0(N)  time. 

Theorem  flft:  Let  S  be  a  set  of  N  points  on  a  circle  with  center  C,  L  be 
a  line  tangent  to  that  circle,  and  S'  be  the  set  of  N  red  and  green  points 
obtained  by  gnomonically  projecting  S  onto  line  L.  There  exists  a 
semicircle  that  contains  all  of  the  points  of  S  iff  the  red  and  green  points 
of  set  S'  are  separable. 

Proof:  Part  (1):  Assume  that  there  exists  a  semicircle  T  that  contains 
the  points  of  S.  Let  the  endpoints  of  semicircle  T  be  called  T1  and  T2 
where  T1  is  on  the  lop  half  of  the  circle  and  T2  is  on  the  bottom  half. 
Since  T1  and  T2  are  diametrically  opposite  each  other,  they  project  to 
the  same  point  T1 '  *  T2'  of  line  L.  We  will  show  that  point  TV  =  T2' 
separates  the  red  and  green  points  of  L.  Any  upper  circular  point  to  the 
left  of  point  T1  maps  to  a  red  point  that  is  to  the  right  of  T1'.  Similarly, 
any  upper  circular  point  to  the  right  of  point  T1  maps  to  a  red  point  that  is 
left  of  TV.  But  for  lower  parts  of  the  circle,  left  and  right  are  not 
reversed.  A  lower  circular  point  to  the  left  of  12  maps  to  a  green  point  to 
the  left  of  T2  and  a  lower  circular  point  to  the  right  of  T2  maps  to  a  green 
point  to  the  right  of  T,?'.  Since  semicircle  T  contains  all  N  points  of  S, 
there  are  only  two  cases  to  consider:  (1)  all  upper  points  of  S  are  to  the 
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left  of  T 1  and  all  lower  point.-,  are  to  the  left  of  T2,  or  (2)  all  upper  points 
are  to  the  right  of  T1  and  all  lower  points  are  to  the  right  of  T2.  For  each 
of  these  two  eases  it  is  easy  to  see  by  the  above  observations  that  the 
point  11'  =  T2‘  .separates  the  red  and  green  points  of  L. 

Part  (2):  Assume  that  the  red  and  green  points  of  L  are  separable.  Let 
point  T*  be  any  point  on  L  that  separates  the  red  and  green  points.  There 
are  two  circular  points,  T1  and  12,  that  project  to  point  T’.  To  complete 
our  proof  we  must  show  that  one  of  the  two  semicircles  determined  by  T1 
and  T2  contains  all  of  the  points  of  S.  We  omit  the  details  because  they 
are  very  similar  to  the  case  analysis  described  in  Port  (1)  above.  □ 

Theorem  Given  N  points  on  a  circle,  we  can  determine  in  0(N)  time 
if  there  exists  a  semicircle  that  contains  all  M  points. 

Proof:  P,y  Theorem  So  we  can  transform  the  problem  of  determining 
inclusion  in  a  semicircle  to  a  problem  of  determining  separability  of  red 
and  green  points  on  a  line.  Since  we  can  determine  if  the  red  and  green 
points  are  separable  by  simply  comparing  the  min  and  max  of  the  red  and 
green  points,  we  can  solve  the  two-dimensional  problem  in  0(N)  time.  □ 


6.2.2.  The  Three-Dimensional  Case 

The  gnomonic  projection  is  a  principal  feature  of  the  0(N  log  N)  time  algorithm  for 
determining  it  N  spherical  points  lie  on  a  hemispherical  cap.  The  definition  of  the 
gnomonic  projection  in  three  dimensions  is  a  straightforward  extension  of  the 
transform  in  two  dimensions; 

Let  P  be  a  point  on  a  sphere  with  center  C  and  let  plane  L  be  tangent 
to  the  sphere.  The  gnomonic  projection  of  spherical  point  P  onto  plane  L 
is  the  point  P*  such  that  points  P,  P',  and  C  are  collinear. 

This  transform  maps  great  circles  on  a  sphere  to  straight  lines  on  a  plane,  allowing 

us  to  use  known  algorithms  for  linear  objects  in  the  Euclidean  plane  to  solve  our 

spherical  problem. 

Wc  assign  the  labels  red  and  green  to  the  projections  on  plane  L  as  in  the 
two-dimensional  case.  If  line  segment  PP'  contains  point  C,  then  we  label  P*  red. 
Otherwise  we  label  P'  green.  Wc  now  use  gnomonic  projection  to  transform  our 
three-dimensional  problem  to  a  two-dimensional  problem  that  we  can  solve  in 
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Figure  G-7:  Gnomonic  projection  of  points  on  a  sphere. 

0(N  log  N)  time. 

Theorem  40:  I  et  S  be  a  set  of  N  points  on  a  sphere  with  center  C,  L  be 
a  plane  that  is  tangent  to  the  sphere,  and  S'  be  the  set  of  red  and  green 
points  obtained  by  gnnmonically  projecting  S  onto  L.  There  exists  a 
hemispherical  cap  that  contains  all  of  the  points  of  S  iff  the  sets  of  red 
and  green  points  of  S'  are  separable  by  a  line. 

Proof:  The  details  are  analogous  to  the  proof  for  the  two-dimensional 
case.  Simply  replace  circle  by  sphere,  semicircle  by  hemisphere,  etc.  □ 


Figure  G-G:  Example  of  nonsepnrahlc  sets  of  red  and  green  points. 
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We  have  just  reduced  the  problem  of  determining  if  a  set  of  spherical  points  lie  in 
a  hemisphere  to  the  problem  of  determining  if  two  planar  sets  of  points  are 
separable  by  a  straight  line.  Shames  and  I  lacy  [94]  present  an  0(N  log  N)  time 
solution  to  the  pi, mar  separability  problem.  They  first  moke  the  observation  that  two 
planar  sets  of  points  are  separable  by  a  line  iff  the  com/ox  hulls  of  the  two  sets  are 
separable  by  a  line.  The  convex  hulls  are  separable  iff  their  intersection  is  empty. 
(See  Figure  G-8  for  an  example  of  nonseparable  sets  of  points.)  Since  we  can 
construct  the  two  convex  hulls  in  G(N  log  N)  time  ([43])  and  then  intersect  them  in 
O(N)  more  time  ([04]),  we  can  determine  separability  of  two  planar  sets  of  0(N) 
points  in  0(N  log  N)  time.  Since  the  gnomonic  projection  of  N  spherical  points  to  N 
planar  points  costs  only  0(N)  time,  we  have 

Theorem  41:  Given  a  set  S  of  N  spherical  points,  we  con  determine  in 
0(N  log  N)  time  whether  or  not  all  N  points  of  S  can  be  fit  into  a 
hemispherical  cap. 

6.2.3.  The  Four-Dimensional  Case 

If  S  is  a  set  of  N  points  on  a  four-dimensional  hypcrsphcrc,  then  in  0(N  log  N)  time 
we  can  determine  if  the  *1  points  of  S  can  be  enclosed  in  a  fGur-dimensional 
hemispherical  cap.  The  algorithm  is  analogous  to  the  algorithm  for  the 
three-dimensional  case.  We  map  the  il  spherical  points  to  N  red  and  green  points  in 
Euclidean  three-space  by  a  (four-dimensional)  gnomonic  projection.  The  N 
hyperspherical  points  fit  on  a  hemi-hypersphere  iff  the  red  and  green  points  are 
separable  by  a  plane.  We  can  determine  separability  of  the  red  and  green  points  in 
0(N  log  N)  time.  The  first  step  is  to  construct  the  convex  hulls  of  the  red  and  green 
points  in  0(N  log  N)  time  ([33]),  and  then  intersect  the  two  convex  hulls  in 
0(N  log  N)  time.  We  can  construct  the  intersection  by  the  algorithm  of  Muller  and 
Preparnta  [73]  or  by  any  0(N  log  N)  lime  algorithm  for  intersecting 
(three-dimensional)  half-spaces.  If  the  intersection  of  the  two  convex  hulls  is 
empty,  then  the  N  hyperspherical  points  of  S  can  be  enclosed  by  a 
hemi-hypcrsphorical  cap. 
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6.3.  Summary 

This  diopter  piesenls  two  dissimilar  sets  of  problems  and  techniques.  For  the 
approximate  Fnr.lidenn  diameter  of  a  planar  set  of  points  we  have  been  interested  in 
a  fast  approximate  convex  hull  algorithm.  One  approach  is  to  successively  use 
rotation  to  approximate  the  fuclidean  metric,  obtaining  an  0(N/(<:1/-))  time  algorithm. 
A  faster  algorithm  can  Lie  obtained  by  a  partitioning  approach;  in  this  case  we  used 
wedge-shaped  "pie-slices"  or  "buckets,"  yielding  an  0(N  +  1/()  time  algorithm.  This 
algorithm  avoids  an  extra  factor  of  0(log  I/O  by  using  the  floor  function  to  drop 
points  into  the  appropriate  "slice." 

To  determine  if  there  exists  a  hemisphere  that  contains  all  N  of  a  given  set  of 
spherical  points  we  found  it  useful  to  apply  a  gnomonic  projection.  This  is  because 
a  hemisphere  is  bounded  by  a  great  circle  and  the  gnomonic  projection  maps  great 
circles  on  a  sphere  to  lines  in  the  Euclidean  plane.  This  enables  us  to  use  a  fast 
planar  algorithm  (separability  of  planar  sets  of  points)  to  solve  the  spherical  problem 
in  0(N  log  N)  time.  Similarly,  N  points  on  a  four-dimensional  sphere  can  be  mapped  to 
Euclidean  three-space  to  obtain  a  problem  (separability  of  two  three-dimensional 
sets  of  points),  which  we  can  solve  in  0(N  log  11)  time. 
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7.  Conclusion 

7.1 .  Transforms  r.nd  Techniques 

This  thesis  presents  a  new  technique  for  construction  of  fast  geometric 
algorithms  --  the  use  of  geometric  transforms.  Although  many  different  kinds  of 
transforms  are  presented,  most  fall  into  one  of  two  classes 

1.  transforms  to  convert  problems  that  are  expressed  in  terms  of  circles 
or  spheres  to  problems  that  are  expressed  in  terms  of  lines  or  planes, 
and 

2.  transforms  to  convert  problems  that  are  expressed  in  terms  of  flats  to 
problems  that  are  expressed  in  terms  of  points. 

Transforms  in  the  first  class  map  points  to  points  and  transforms  in  the  second  class 
are  duality  transforms.  In  Appendix  III  we  list  the  transforms  used  in  this  thesis, 
their  important  properties,  and  some  of  their  applications.  In  addition,  Table  7-1 
summarizes  "typical  uses"  for  four  of  the  transforms  that  have  been  particularly 
useful. 
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Projection  (orthographic  or  "radial") 

Solve  a  K-dimensional  problem  ns  n  (K-1  )-dimensional  problem. 
(Reducing  dimensionality  typically  simplifies  a  problem.) 

Embedding  Solve  a  K-dimensional  problem  as  a  (K+1  )-dimensional  problem. 

(The  extra  degree  of  freedom  can  sometimes  allow  an 
interpretation  of  the  problem  which  is  not  possible  in  K-space.) 

Point  /  flat  Duality  Solve  a  problem  involving  flats  as  a  problem  involving  points.  If 
half-spaces  (rather  than  flats)  are  involved,  then  the 
t  half-spacer,  naturally  partition  into  two  sets  (UPPER  and 

I  OWI-ft)  since  each  flat  determines  two  half-spaces. 

Inversion  (Stcreogrnphic  Projection) 

Converts  problems  involving  circles  (spheres)  to  problems 
involving  lines  (planes).  Since  circles  (spheres)  are  intimately 
related  to  the  Euclidean  metric,  inversion  may  be  useful  in 
problems  involving  the  Euclidean  distance  between  points. 

Table  7-1:  Heuristics  for  use  of  some  transforms. 


7.2.  New  Results 

In  the  preceding  chapters  we  hove  presented  several  examples  of  the  use  of 
geometric  transforms.  Some  of  the  major  results  are: 

-  We  relate  the  Euclidean  Voronoi  diagram  of  a  set  of  points  to  the 
convex  hull  of  a  transformed  set  of  points,  yielding  an  0(N  log  N)  time 
planar  algorithm,  which  is  optimal  to  within  a  constant  factor. 

-  We  transform  the  problem  of  finding  planes  of  support  of  a  convex 
polyhedron  to  the  problem  of  finding  all  overlapping  regions  among  two 
straight-line  planar  graphs.  This  gives  an  0((N  +  K)  log  N)  time  and  0(N) 
storage  algorithm  for  computing  the  Euclidean  diameter  of  a  set  of  N 
points  in  three-space,  where  K  is  the  number  of  pairs  of  antipodal 
vertices  nf  the  convex  hull  of  the  N  points. 

-  We  relate  the  intersection  of  UPPER  half-spaces  to  the  convex  hull  of 
a  finite  set  of  points.  This  leads  to  an  0(N  log  N)  lime  algorithm  for 
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constructing  the  intersection  of  N  UPPER  three-dimensional 
half-spaces,  which  is  optimal  to  within  a  constant  factor. 20  We  apply 
the  floor  function  to  obtain  an  0(M+1/O  time  algorithm  for 
approximating  the  diameter  of  N  planar  points  to  within  a  factor  of 
1  +  (.  (This  and  Bentley,  Faust,  and  Preparala's  [7]  0(N  +  I/O  time 
planar  algorithm  are  the  best  known  algorithms  for  approximating  the 
Cuclidean  diameter.) 

-  We  transform  the  union  of  N  planar  disks  to  an  intersection  of  N 
three-dimensional  half-spaces,  yielding  an  0(N  log  N)  time  algorithm  for 
constructing  the  union  of  N  planar  disks. 

-  There  are  also  a  number  of  minor  hut  new  results: 

*  Location  of  an  arbitrary  point  in  a  set  of  N  planar  disks  in 
0(log  N)  lime  with  only  0(N  log  M)  preprocessing  time  and  0(N) 
storage  (Section  3.2.5). 

*  0(.M  log  N)  time  algorithms  for  constructing  closest  and  farthest 
point  Voronoi  diagrams  for  N  spherical  points  and  for  the  sides  of 
a  convex  N-gon  (Sections  4.2  and  4.3). 

*  Determining  in  0(N  log  N)  time  whether  N  points  on  a 
four-dimensional  sphere  can  be  fit  in  a  hemisphere  (Section 
6.2.3). 


7.3.  Open  Problems 


There  are  still  several  problems  that  remain  unsolved.  The  list  below  describes 
several  of  the  major  outstanding  problems. 

-  One  of  the  major  unsolved  problems  of  this  thesis  is  to  prove  an 
Q(N  log  N)  lime  lower  bound  for  computing  the  Euclidean  diameter  of  N 
planar  points  under  a  model  of  computation  strong  enough  to  compute  it 
in  0(N  log  N)  lime.  (Appendix  II  describes  on  approach  toward  such  a 
proof.)  Shamos  [02]  conjectures  that  for  any  metric  with  a 
continuously-turning  tangent  the  worst-case  lower  bound  is  Q(N  log  N) 
time. 


result,  nlltnwjh  rimil;,r  to  the  work,  of  Pre;.iVata  and  Mjllcr  [o-j,  was  done  independently. 
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-  What  is  the  expected  number  of  pairs  of  antipodal  vertices  among  N 
three-dimensional  points  (under  any  interesting  distribution  of  points)? 
(Thin  determines  the  expected-lime  of  the  Euclidean  diameter 
algorithm  of  Section  f>.?.) 

-  For  an  approximate  diameter  vie.  and  Dentiay.  Faust,  and  Preparata 
[7]  have  shown  an  0(N  +  I/O  time  upper  bound,  but  the  lower  bound 
remains  an  open  question. 

-  The  "bucket  transform"  of  Weidc's  0(N)  expected-timc  sort  [99]  has 
already  been  extended  to  0(N)  expcolcd-time  planar  Voronoi  diagrams 
by  Bentley,  Weide,  and  Yao  [IS].  The  power  of  this  transform  derives 
from  the  use  of  the  floor  function.  Most  algorithms  and  proofs  for 
lower  bounds  use  models  of  computation  that  involve  only  comparisons 
between  analytic  functions  of  the  input,  but  these  results  are 
becoming  dated  as  the  floor  function  is  applied  to  more  problems. 
Some  other  notable  uses  of  the  floor  function  are  the  linear-time 
largest  gap  algorithm  of  Gonzalez  [47],  linear  expected-time 
closesl-poinl  algorithms  of  Yuval  [104]  and  Rabin  [S6],  and 
0(N  log  log  M)  worst-case  time  closest-poir  algorithm  of  Fortune  and 
llopcrofl  [41].  These  algorithms  partition  their  problem  into  "buckets" 
and  use  the  floor  function  to  compute  quickly  for  each  point  the  bucket 
into  which  it  should  he  dropped.  What  other  uses  can  we  find  for  the 
floor  function  in  geometric  algorithms? 

-  The  algorithms  for  intersection  of  Imlf-spaccs,  planar  Euclidean  Voronoi 
diagrams,  spherical  Voronoi  diagrams,  and  diameter  of  a  set  of  points 
all  use  Iho  convex  hull  of  a  set  of  points.  (Also,  the  algorithms  that 
involve  an  intersection  of  half-spaces  indirectly  use  convex  hulls  of 
sets  of  points.)  Shnmos  [91]  describes  how  a  convex  hull  is  used  in 
isotonic  regression  and  Silverman  and  Titterington  [96]  use  a  convex 
hull  to  find  the  smallest  covering  ellipse  of  a  planar  set  of  points. 
What  other  uses  can  be  found  for  convex  hulls? 


7.4.  Conclusion 

Transforms  are  often  useful  for  converting  apparently  difficult  problems  to 
instances  of  problems  that  are  solvable  by  well-known  methods.  !n  this  thesis  we 
have  presented  a  set  of  techniques  for  applying  geometric  transforms  to  geometric 
problems  that  provides  another  sel  of  tools  for  the  designer  of  fast  geometric 


1 


24  Decernhci  1  070. 


Geometric  Transforms 


PAGE  123 


algorithms.  Curllx  rmore,  in  tin-  prc.ouss  of  illustrating  the  application  of  these 
transforms  to  sovor.il  problems,  wo  have  obtained  several  useful  and  interesting 
algorithms. 
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Appendix  I 

Finding  a  Good  Orientation  for  Flats 

One  of  the  major  problems  v/itli  representing  lines,  planes,  or,  in  general,  flats  In 
slope-intercept  form  is  that  vertical  flats  cannot  be  represented  and  near-vertical 
flats  cause  large  rouiul-off  error.  If  vertical  and  near-vertical  flats  can  be  avoided, 
however,  the  slope-intercept  form  is  very  convenient  because  of  the  properties  of 
the  transform 


K-1 

xs  =  2  'V,  +  aK  (a„  a2 . aK) 

i=1 

as  described  in  Section  3.1.  In  this  appendix  we.  describe  how  rotation  of  the  flats 
can  help  avoid  tin:  occurrence  of  vertical  or  near-vertical  flats. 

There  ore  two  cases  to  consider: 

(1)  The  restricted  case:  If  the  relation  is  being  performed  to  enable  an 
intersection  of  N  UPPER  half-spaces  (as  in  Section  3.1),  then  we  must 
ensure  that  after  rotation  all  N  half-spaces  remain  UPPER  half-spaces 
(rather  than  LOWER  half-spaces). 

(2)  The  general  case:  In  general  there  is  no  restriction  as  in  Case  (1). 

1.1.  Case  (1) 

Since  we  want  to  determine  what  engles  to  rotate  the  flats  through,  it  will  be 
useful  first  to  represent  the  flats  by  their  angles  rather  than  slopes.  Specifically, 
we  map  the  slope-intercept  representation  to  an  angle-of- inclination  representation 
by 


K-1 

xK  =  X  tan(0,)x,  +  aK  -*  (0V  02 . 0K.t) 

i=1 


where  0j  *  [-n/2,  n/2),  j  *  1 . K-1.  In  the  following  discussion  a  K-dimensional 
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problem  of  flats  will  therefore  be  treated  ns  a  K-1 -dimensional  problem  of  points.  A 
flat  will  be  called  near-vertical  if  one  of  the  K  -  1  corresponding  angles  0j  is  within 
a  given  small  angle  <-  of  being  vertical  (n  / 2  or  -Tt/2).  We  will  also  use  the  notation 
for  the  jth  coordinate  (angle)  for  the  ith  flat. 

We  want  to  find  a  rotation  (<yj,  a 2 . sucI> 

-TT/2  +  <  <  ffj  +  Qjj  i  rt/2-f.  for  i  =  1,  . . .,  N.  j  s  1 . K-1.  (33) 

We  can  easily  find  such  a  rotation,  if  it  exists,  in  0(KN)  time  and  storage.  First 
compute  the  maximum  and  minimum  angles  for  all  K  -  1  coordinates 

max  ,  min .  .  .  . 

maxj  »  .  ©ij  and  min-  =  0^,  for  j  *  1, . . .,  K-1, 

In  O(KN)  time  and  then  determine  if 

|  It  -  (maxj  -  minj)  |  >  2(  for  j  *  1 . K-1.  (34) 

If  Condition  (3d)  is  satisfied  then  a  rotation  satisfying  Condition  (33)  exists.  One 
such  rotation  j ,  . . ..  )  is 

*  -  (maxj  ♦  minp  /  2,  j*1 . K-1. 

It  may  not  be  necessary  to  do  any  rotation  at  all;  if 

|  tt/2  -  maxj  |  >  *  and  j  IT/2  +  roinj  |  2  €  for  j  *  1, . . .,  K-1  (35) 

then  no  rotation  is  required.  Tiie  probability  that  Condition  (35)  is  satisfied  depends 
on  the  probability  distribution  of  the  points  (Ojj,  0 ^  •  • ..  Gj.K-i^*  the  density  of 
the  distribution  is  uniform  (and  ell  points  are  independent)  tiien  Condition  (35)  will  be 
satisfied  with  probability 

•  ■  klfV-1  \ 

P(  No  rotation  rcqtiircd  )  =  (-1— 

\  If 

Condition  (34)  (existence  of  a  rotation  satisfying  Condition  (33))  is  satisfied  with 
probability 

P(rotation  for  Case  (1)  exists)  *  P(maXj-minj  i  Tt-20,  j  *  1, . .  .,K-1  ■ 

«  P(max^-min1  S  IT-20*’1 
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wnere 


P(max,-min,  S  n-20  =  f*  f*2 

'  ,  J  X2*0  J  x  i  *i 


r*2 


a2 


max(0,.»2-r+^i)  dx,dx2 


-FPT 


dx,  dx2 


/IT-2*1  N  2(N 

\  n  J  it  1  rt  j  ’ 


1.2.  Coco  (2) 

If  we  allow  arbitrary  rotations  then  the  K-1  -dimensional  points  (01te2 . Gk-i) 

may  "wrap  around"  the  boundaries  of  the  cube  of  angles  [-TT/2,Tt/2)^~\  Whereas 
the  analysis  of  Case  (1)  required  only  computations  related  to  the  max's  and  min's 
for  each  of  the  K  -  1  coordinates  (independently),  the  analysis  of  Case  (2)  involves 
the  "largest  great-circle-shaped  gap"  between  N  points  on  a  K-dimensional  sphere. 

The  mapping  that  we  use  is  (almost)  the  point  /  flat  duality  of  Zolnowsky  [106].  A 
flat  is  first  translated  so  that  it  is  tangent  to  the  unit  sphere  (centered  at  the 
origin)  anti  the  transform  is  the  point  of  tangency.  Since  there  are  two 
(diametrically  opposite)  points  at  which  a  fiat  of  a  given  orientation  can  be  tangent, 
there  are  two  possible  points  to  which  a  flat  can  map.21  We  can,  of  course, 
describe  this  transform  in  cartesian  coordinates,  but  when  we  wish  to  speak  of  the 
angles  of  inclination  it  will  be  more  useful  to  use  the  mapping 

K-1 

XK  s  Z  •cotO.JXj  ♦  aK  -»  (0,,e2 . e,..,).  (36) 

1*1 

4  * 

where  the  points  •  •  •»  0K-l)  ar®  interpreted  to  be  normal  geodesic 

coordinates  of  a  K-1 -dimensional  sphere.  (Each  K-1 -tuple  represents  two 


21 

Wo  can  also  c>  press  tlx?  transform  as  the  point  /  flat  duality  of  Section  3.1  followed  by  orthographic 
projection  to  the  planes  *  *1,  reversing  Uie  sign  of  all  coordinates,  and  tl\en  gnomonie  projection  to  the  surface 
of  tho  sphere. 
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diametrically  opposite  points  on  a  sphere.)  In  each  coordinate  system  this  transform 
maps  each  flat  to  a  point  on  the  unit  sphere  such  that  the  vertical  flats  map  to 
points  on  the  groat  circle  0,  =  ff/2,  i  =  1,....  K-1  (or,  in  cartesian  coordinates, 
where  the  sphere  intersects  the  flat  x^  =  0).  The  near-vertical  flats  map  to  points 
near  this  great  circle.  We  will  present  an  analysis  for  two  and  three  dimensions. 

1.2.1.  Case  (2)  for  Lines  in  a  Plane 

Lines  in  the  plane  transform  by  Equation  (3G)  to  points  on  the  top  half  of  the 
circle  *?-  +  y^  =  1  and  their  diametrical  opposites  on  the  bottom  half.  A  line  L  is 
vertical  iff  it  maps  to  the  points  (x,y)  =  (1,0)  (Q  =  0)  and  (x.y)  =  (-1,0)  0  =  L  is 
near-vertical  iff  it  maps  to  points  within  an  angle  <r  of  (1,0)  and  (-1,0).  There  exists 
a  rotation  of  a  set  of  lines  S  such  that  none  are  near-vertical  iff  there  exists  a 
rotation  of  their  corresponding  points  S'  on  the  circle  such  that  none  are  within  an 
angle  (  or  (1,0)  or  (-1,0).  This  occurs  iff  there  exists  a  rotation  of  the  points  (1,0) 
and  (-1,0)  such  that  lliey  are  not  within  an  angle  £  of  any  of  the  N  points  of  S’, 
which  occurs  iff  there  exists  a  gap  of  at  least  Z<  between  adjacent  points  of  S'  on 
the  circle.  We  can  determine  this  in  0(N)  time  by  (a  modification  of)  Gonzalez's 
largest  gap  (algorithm  [47].  (It  should  be  noted,  however,  that  Gonzalez's  algorithm 
uses  the  floor  function.  If  only  analytic  functions  are  used,  then  the  best  known 
algorithm  for  the  largest  gap  problem  takes  0(N  log  N)  time.) 

1.2.2.  Case  (2)  Tor  Planes  in  Throe-Space 

Planes  in  three-space  transform  by  Equation  (3G)  to  points  on  the  top  half  of  the 
sphere  +  y^  +  z^  =  1  and  their  diametrical  opposites  on  the  bottom  half.  Vertical 
planes  map  to  points  on  the  great  circle  0is02  =  W2  (where  the  sphere 
intersects  the  plane  ?  -  0)  and  near-vertical  planes  lie  within  an  angle  <  of  this 
great  circle.  There  exists  a  rotation  for  a  set  of  planes  S  such  that  no  plane  is 
near-vcrticol  iff  there  exists  a  rotation  for  this  great  circle  such  that  no  point  is 
within  an  angle  <  of  any  of  the  points  of  S'.  Whereas  in  two  dimensions  we  obtained 
a  problem  that  involved  the  (angular)  distance  between  pairs  of  points  on  a  circle,  in 
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three  dimensions  we  hove  a  problem  that  involves  the  (angular)  distance  between 
points  and  a  great  circle  on  a  sphere. 

Searching  for  a  great  circle  that  satisfies  our  conditions  seems  more  difficult  than 
searching  for  a  point  because  a  great  circle  is  a  nonlocal  object.  We  can,  however, 
apply  a  point  /  great  circle  duality  to  obtain  an  equivalent  point-searching  problem. 
The  transform  is  simple; 

spherical  point  -»  great  circle  farthest  from  the  point,  and 
great  circle  -*  the  two  points  farthest  from  the  great  circle. 

Since  this  duality  preserves  the  angular  distance  between  a  point  and  a  great 

circle,  a  spherical  point  P  and  a  great  circle  C  are  an  angular  distance  of  £  apart  iff 

their  transforms  are  *  opart.  Our  new  problem  is  thus 

Given  N  great  circles  on  a  sphere,  find  a  point  on  the  sphere  that  is  not 
within  an  angle  i  of  any  of  the  great  circles. 

This  problem  can  he  solved  by  the  spherical  analog  of  a  planar  nearest  line  Voronoi 

diagram.  Unfortunately,  since  the  N  great  circles  partition  the  sphere  into  Q(N^) 

regions  this  approach  will  take  at  least  Q(fj2)  time.  In  the  worst  case,  it  can  take 

longer  to  find  a  good  orientation  for  a  set  of  N  planes  than  it  takes  to  intersect  f4 

three-dimensional  half-spaces.  If  f  is  small  compared  to  N,  however,  we  might  be 

able  to  find  a  good  orientation  in  0(N)  expected-time. 
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Appendix  !! 

Relation  of  Diameter  to  Empty  Intersection 


In  this  appendix  we  describe  a  relationship  between  the  Euclidean  diameter  and 
on  empty  intersection  problem  that  may  lead  to  an  Q(N  log  N)  time  lower  bound  for 
the  diameter  problem.  The  empty  intersection  problem  that  we  are  interested  in  is  a 
variation  of  the  problem  of  Reingold  [S7]: 

Given  two  sets,  P  and  0,  each  of  N  reals  in  the  interval  [0,1), 
determine  if  there  exist  p(  P  and  q  (  0  such  that  p  =  q. 

We  convert  the  empty  intersection  problem  into  a  diameter  problem  by  using 

inversion  (Section  3.2.3).  The  inversion  transform  is  determined  by  two  parameters, 

the  center  of  inversion  and  the  radius  of  inversion.  Wc  shall  choose  the  origin  as 

the  center  of  inversion  and  let  the  radius  of  inversion  be  one.  The  transform  is  thus 


(x.y)  -» 


’  x 

\x  +  y 


y- 


For  our  purposes  the  only  important  property  of  this  transform  is  that  it  transforms  a 
line  that  does  not  pass  through  the  origin  into  a  circle  that  does  pass  through  the 
origin.  In  particular,  the  line  y  =  1  maps  to  a  circle  with  radius  1  /2  centered  at  the 
point  (0,1/2). 

0  1 


Figure  7-1:  Mapping  the  points  of  sets  P  and  Q  to  a  circle. 
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Our  construction  first  mops  the  N  elements  of  set  P  to  points  on  the  line  y  =  1  and 
then  applies  inversion  to  map  those  points  to  points  on  the  circle 
X2  +  (y  -  1/2)2  =  (1/2)2: 


P 


(P.1) 


o  *  *"> 

(1  +  P‘  1  +  P‘, 


p  6  P. 


Note  that  since  the  elements  of  P  are  chosen  from  the  set  [0,1),  all  of  the  N  points 
generated  from  P  will  lie  within  an  arc  of  only  one  quarter  of  the  circumference  of 
the  circle.  (See  Figure  7-1.)  Our  transform  for  the  elements  of  set  Q,  on  the  other 
hand,  generates  points  on  a  quarter-circle  arc  that  is  diametrically  opposite  the 
points  from  set  P.  More  precisely,  since  we  want  to  convert  the  empty  intersection 
problem  to  a  diameter  problem  we  choose  the  transform  for  the  elements  of  set  Q  so 
that  identical  elements  of  sets  P  and  Q  will  be  diametrically  opposite  on  the  circle. 
The  transform  for  set  0  is  thus 


<1 


(-q.  0) 


-q  q2  \ 

1  +  q2  *  1  +  q2j  ’ 


q  €  Q. 


Since  all  of  the  elements  of  sets  P  and  0  are  transformed  to  points  on  a  circle,  the 
maximum  possible  diameter  of  the  2N  points  is  the  diameter  of  the  circle.  This 
maximum  is  achieved  iff  two  of  the  2M  points  are  diametrically  opposite,  which 
occurs  iff  there  exists  a  p  (  P  and  q  <  0  such  that  p  =  q.  We  have  just  proved 

Theorem  d2:  Given  a  model  of  computation  strong  enough  to  support 
the  inversion  transform,  we  can  transform  an  instance  of  the  empty 
intersection  problem  into  an  N-point  diameter  problem  in  0(N)  time. 

To  apply  our  lower  bound  for  the  empty  intersection  problem  to  the  diameter 
problem,  wc  must  prove  the  lower  bound  under  a  model  of  computation  strong  enough 
to  support  our  construction  for  solving  the  empty  intersection  problem  as  a  diameter 
problem.  Roingold's  Q(N  log  N)  lower  bound  unfortunately  allows  only  comparisons 
with  linear  functions  of  the  input,  which  is  not  strong  enough  for  our  purposes.  One 
model  of  compulation  that  is  strong  enough  is  a  decision  tree  that  allows  only 
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comparisons  between  multivariate  polynomials  of  degree  at  most  K  (for  some  fixed  K 
>  4)  at  internal  nodes.  It  is  important  that  the  degree  be  bounded  because 
otherwise  we  could  solve  the  empty  intersection  problem  in  constant  time  by  making 
the  comparisons 


N  N 

n  (p.  ■  <ij)  -  ° anci  n  (Pi  *  Hj)  *  o. 

i.j  i.j 

Comparisons  between  bounded  degree  polynomial  functions  are  sufficiently  strong 
that  we  can  support  our  construction  for  solving  empty  intersection  as  a  diameter 
problem. (Furthermore,  we  can  construct  0(M  log  N)  time  diameter  algorithms 
under  this  model  of  compulation.)  It  remains  an  open  problem  whether  or  not  we  can 
prove  an  ft(N  log  N)  time  lower  bound  for  the  empty  intersection  problem  under  this 
model  of  computation. 


22lt  may  at  first  appear  that  wo  cannot  compute  the  inversion  transform  of  our  construction  with  finite  degree 
polynomial  functions  because  it  involves  division.  That  is,  in  a  strict  sense,  true,  but  we  can  simulate  it  by 
representing  a  quotient  p/q  as  an  ordered  pair  (p,q). 
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Appendix  i!i 

Georneiric  Transforms  and  Applications 

In  this  appendix  we  summarise  the  transforms  used  in  this  thesis  (or  used  for 
geometric  problems  not  described  in  this  thesis)  and  list  properties  of  and 
applications  for  each.  The  transforms  are  collected  into  three  categories: 
Point-to-Point  Transforms,  Duality  Transforms,  and  Miscellaneous  Transforms. 

111.1 .  Point-to-Point  Transforms 

Most  of  the  point-to-point  transforms  fall  into  one  of  tivo  classes;  continuous 
and  invertible  transforms  and  projections.  All  continuous  and  invertible 
transforms  f  arc  potentially  applicable  to  union  or  intersection  problems 
because  they  satisfy 

f(AV  D)  -  f(A)U  f(D)  and  f(AC\B)  =  f(A)fU(B). 


Rotation 

Preserves  size  and  shape  while  altering  orientation. 

-  Convert  1.0J  diameter  to  L1  diameter  in  the  Euclidean  plane  (Section 

2.6). 

General  linear  transform 

Maps  linear  quantities  to  linear  quantities  while  stretching  and  rotating. 

-  Generalizes  problems  with  reclilincarly  oriented  lines  or  line  segments 
to  problems  with  lines  or  line  segments  at  two  arbitrary  angles 
[15,  12]. 

-  Ono-dimensional  Johncon-Mchl  crystal  growth  model  [4G]  transforms  to 
two-dimensional  maxima  of  vectors. 

-  Derivation  of  point  /  flat  duality  from  inversion  and  a  linear  transform 
(Section  3.3). 


Orthographic  projection 
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Hcdu cos  dimensionality  by  eliminating  one  of  the  cartesian  coordinates  and 
leaving  the  others  unaffected. 

-  Searching  a  Voronoi  diagram  [89,  33,  31]. 

-  Multi-dimensional  dividc-and-conquer  algorithms  that  solve  lire  merge 
step  of  a  K-dimensional  problem  as  a  K-1 -dimensional  problem  [3]. 

-  Nearest  (farthest)  side  diagram  of  a  convex  polygon  (Section  4.3). 

-  Euclidean  diameter  of  N  points  in  three-space  in  0((N  +  K)  log  N)  time, 
where  K  is  the  number  of  pairs  of  antipodal  vertices  on  the  convex  hull 
(Section  5.2). 

-  Three-dimensional  Chcbyshev  regression  (Section  5.2.3). 

-  Least-squares  isotonic  regression  [91]. 

-  Least  squares  regression. 

Perspective  transformation 

Lines  map  to  lines  and  pianos  map  to  planes  while  preserving  perspective 
information.  Maps  K-dimcnsionat  rays  to  K~1 -dimensional  line  segments. 

-  Transforms  a  perspective  projection  (visibility)  problem  to  an 
orthographic  projection  problem  [97]. 

Radial  projection  (about  a  point) 

Preserves  spherical  angle  while  reducing  dimensionality. 

-  Dual  of  spherical  Voronoi  diagram  (from  convex  hull)  (Section  4.2). 

-  Spherical  Voronoi  diagram  (from  intersection  of  half-spaces)  [22]. 

Gnomonic  projection 


A  two-lo-one  transform  that  maps  great  circles  on  a  sphere  to  lines  in  the 
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-  Determine  if  N  spherical  points  can  fit  in  a  hemispherical  cap  [105] 
(Section  G.2). 

Stereographic  projection 

A  conformal  mapping,  bul  il  reducer,  the  sense  of  the  angle . 

-  This  is  a  special  ease  of  inversion  in  three-space.  (See  inversion.) 

Inversion 

i 

Inversion  is  a  circular  transform;  [circles  map  to  circles  (where  a  fine  is 
considered  to  be  a  circle  of  infinite  radius).  In  particular, 

-  A  circle  that  passes  tluougli  the  center  of  inversion  maps  to  a 
line  that  does  not  pass  through  the  center  of  inversion. 

-  The  interior  of  a  circle  that  contains  the  center  of  inversion 
maps  to  the  exterior  of  a  circle  that  contains  the  center  of 
inversion. 

-  Other  properties  are  described  in  [36], 

-  Union,  intersection,  subtraction  of  disks  (Section  3.2). 

-  Nearest  (farthest)  point  Voronoi  diagram  (Section  4.1)  [23]. 

-  Nearest  (farthest)  edge  diagram  (Section  4.3). 

-  Derivation  of  point  /  flat  duality  from  inversion  and  linear  transform 
(Section  3.3). 

-  Mapping  of  points  on  a  line  to  points  on  a  circle  in  Appendix  II. 

"Bucket"  transform  (floor  function) 

A  discontinuous  function  but  available  as  a  primitive  on  most  machines.  Usoful 
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-  O(N)  expected-time  sort  [00], 

-  0{N)  cxpocted-time  convex  hull,  Voronoi  din rjr r.m  [1C]. 

-  0(N  +  I/O  time  (-approximation  for  Euclidean  diameter  of  points 
(Section  G.1 )  [7]. 

Embedding  into  a  higher  dimension 

Adds  a  dogma  of  f random  to  the  problem.  Linos  con  become  planes,  circles 
can  become  spheres,  etc. 

-  Union,  intersection,  subtraction  of  disks  (Section  3.2). 

-  Nearest  (farthest)  point  Euclidean  Voronoi  diagram  (Section  4.1). 

-  Nearest  (farthest)  side  diagram  of  a  convex  polygon  (Section  4.3). 

-  Lower  bounds  for  union  (intersection)  of  disks  or  half-planes  (Section 
3.2)  or  convex  hull  of  points  or  triangulation  of  a  set  of  points 
[91].  (Sort  N  reals  by  mapping  them  to  N  planar  points,  N  disks,  or  N 
half-planes.) 

-  0(N  log  M)  time  algorithm  for  least-squares  isotonic  regression  [01]. 

111.2.  Duality  Transforms 

Useful  for  transforming  problems  involving  nonlocal  objects  (such  as  flats)  to 
problems  involving  points. 

N  Points  In  1  Space  /  1  Point  in  N  Space  Duality 

Transforms  an  N  point  problem  into  a  one  point  problem. 

-  Element-uniqueness  lower  bound  [34]. 

-  Epsilon-closeness  lower  bound  [42]. 

-  Lower  bounds  for  sorting,  insertion,  finding  max  with  analytic  functions 
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[43]. 

Zolnowsky's  plane  /  point  on  unit  sphere  duality 

-  0(N  log  N)  time  intersection  of  H  half-spaces  [10G]. 

-  Finding  a  good  orientation  for  flats  {Case  (2)  of  Appendix  I.) 

Plane  /  Homogeneous  Pluckcrian  coordinate  duality 

A  more  homogeneous  representation  of  the  point  /  flat  duality  below. 

-  0(N  log  N)  time  intersection  of  ft  half -spaces  [84]. 

Slope  Intercept  form  of  Point  /  Flat  duality 

Preserves  x ^  coordinate  distance  between  a  point  and  a  flat,  thereby 
preserving  incidence  between  points  and  flats.  Preserves  above /belowness 
between  points  and  flats. 

-  0(N  lor*  N)  time  intersection  of  N  UPPER  (or  N  LOWER)  half-spaces 
(Section  3.1 ). 

*  Closest  (farthest)  side  of  a  convex  polygon  (Section  4.3). 

*  N  points  on  a  hemisphere:  intersection  of  half-spaces  solution 
(Section  G.2). 

-  Linear  programming  [23]. 

-  Impossible  three-dimensional  scenes  [03], 

-  Diameter  of  N  points  in  thrce-f.pncc  in  0((N  +  K)  log  N)  time  (Section 
5.2). 

-  Least-squares  isotonic  regression  [91]. 

-  Three-dimensional  Chebyshev  regression  (Section  5.2.3). 
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-  Tour  two-dimensional  line  problems:  center  of  lines,  intersection  radius 
of  lines,  minimal  spanning  segment  of  lines,  diameter  of  lines  [21]. 

-  Number  of  exterior  points  in  intersection  of  N  lines  is  linear  [21]. 

Generalized  (Slope-Intercept)  Duality  Transform 

-  Intersection  of  I  lalf-K-Spaces  (Section  3.1.5). 

Circle  /  Point  Duality 

-  Relation  between  inversion,  Point  /  flat  duality,  and  convex  hulls 
(Section  3.3). 

Spherical  Point  (Pair)  /  Great  Circle  Duality 

Preserves  sphcrie.nl  angle  between  a  point  end  a  great  circle. 

-  0(N  log  N)  lime  algorithm  for  spherical  farthest  line  Voronoi  diagram 

[22]. 


III. 3.  Miscellaneous  Transforms 

Point  to  locus  of  a  set  of  points  transforms 

-  Digitization  problem.  Given  M  digitized  points,  find  the  set  of  all  lines 
that  have  that  digitization  [22]. 

-  Inclusion  of  N  points  on  a  hemisphere  --  intersection  of  half-spaces 
solution:  each  point  transforms  to  the  set  of  points  from  which  it  is 
visible  (in  this  case  a  half-space).  (Sec  Section  4.2.) 


-  Transform  polyhedral  obstacles  to  locus  of  forbidden  positions  of  a 
reference  point  of  the  moving  object  [71]. 
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Polygon  to  string  transform 

-  O(N)  time  algorithm  for  similarity  of  polygons  [72]. 
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Dunham,  Kelly,  and  Tolle  81 
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Elemcnt-Uniqnene.  :■  f.,  I". t 

Cnil)»<liiiiiri  into  a  liifiiioi  dimension  Cl,  13.3 

Union  of  cl'  sk  r.r. 

Euclidean  Vorfiii.ii  diagram  09 
Spl.etii.il  Vornnoi  tli,n:r am  7d 
Nearest  0(10"  diagiam  70 

Flat  13 

Floor  function  111,  122,  ijg 
Floyd  C,  73 

Foituno  and  Hop, , oft  S,  17,  122 
Frcdman  and  Woide  13.3 
Friedman  77,  54,  13S 

Carey,  fiiahani,  and  John 'tin  1 1 
General  linear  transform  135 
Gilbert  03.  135 
Gnomcme  projection  110,,  130 
Gorcale,'  17.  113.  133,  123 
Graham  5,  35.  A*,,  33,  105,  m,  117 
Grunbaum  0,  101 

Harary  07 
Hartignn  03,  P3 
Hocking  and  Young  £4 
Hoey  7,  OG 
Huffman  39,  130 

Intersection  of  Dn-ks  52 
see  Union  of  Disks 
Intersection  of  Half-Spaces  23 
representation  2d 

Lower  [Join  id  25 

Union  of  disks  5ft 
Nearest  edge  rli.icr.iai  73 
Choosing  orient-alloy  125 
Intersection  Problems  7,  23 
Inversion  5.1,  C2,  137 
Union  of  disks  55 
derivation  of  point  /  flat  duality  £1 
Voronoi  di.vn  am  09 
Nearest  edge  diagram  70 
Isotonic  regression  122 

Jarvis  0 
Jcroslow  pi 

Johnson  and  Pi  epai  ala  1 12 

Kan.ade  39 
Kelly  81 
Kirkpatrick  63,  75 
Klee  and  Mmty  at 
Kleinrook  12 
Kmilh  18 
Kruskal  9 

Knng,  Luceio,  and  Prepar.yt.y  10 
Lee  76 

Lee  .and  Prep.ar.ata  9.  63,  95 
Lee  and  Wong  63 
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Linear  Programming  31 
Liplon  n  riri  Ini  jam  10,  68,  Co,  95 
Lower  Bound  13 
Convex  Hull  G 

Element  Uniqueness  Problem  3 
Closest  Pmr  3 

M.avimn  of  n  Sol  of  Vectors  10 
Inter  section  of  linlf-Spae.es  30 
Sorting  27 
Union  of  0‘  '  s  53 
Voronoi  Diagram  G3 
Eucbdcan  Diameter  00,  131 
Empty  Inter  section  PioWem  131 
Locano-Pero;  ami  Wesley  130 

Mannclter  1'11 
Max  min  of  vector  $  10 

Medial  Ay. is  of  a  Convex  Polygon  76 
see  Mearest  Edge  Diagram 
Metrics  15 

Mirnni'im  Spanning  Tree  9 
Muller  and  Preparata  1 17 

Neare  d  Edge  Diagram  76 
see,  Voronoi  Diagram 

Orthographic  projection  1.35 
Nearest  ccge  diagram  73 
Linear  programming  32 
Euclidean  diameter  87 

Papadimitnou  and  Steiyliu  1 1 
Parker  12 

Perspective  transformation  13G 
Pol.l  10,  83 

Point  /  Flat  Duality  32,  45.  61.  139 
Intersection  of  Half-Spaces  33 
derivation  53 
Linear  Programming  82 
Euclidean  diameter  37 
Preparata  6,  10,  63,  76,  96 
Preparata  and  Hong  G,  40,  GO,  74,  8°,  117 
Preparata  and  Muller  7,  32,  121,  139 
Prim  9 


nahln  8,  17,  122 
Radial  projection  13G 
Reingold  131 
Rotation  130 
L  j  Diameter  20 

Saxe  and  Bentley  1 1 
Stiamo  and  Heey  52 

Shames  5,  9,  12,  16,  23,  75,  S3,  3-1,  30,  100,  105,  111,  121,  122,  136,  13S,  139 
Sliomor.  and  Hoey  7,  3,  27,  63,  67,  82,  00,  1 17 
Slmmos  and  Vuval  101,  103 
Silverman  arid  Titter ingtou  122 

Spherical  piohlems  (transformed  to  linear  problems)  119 
Union  of  disks  62 
Euclidean  Voronoi  diagram  63 
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NCilfC't  tll-Vf  dill  /  0 

Fitting  points  on  «»  I  iohm  -ptivic  1  12 
Splici  K  Al  VoioiiOi  di.vji  «ii n  /  3 
StCfC'Hiitiphii;  projection  137 
SuUkmIAIhI,  fyioult,  SU‘ HIMeklM  130 

Union  of  C>rV."'  62,  132 

rcprc:»?*i»lnliofi 
Lower  Round  S3 

Van  l**enweri  and  Wood  7 
Voronoi  Dmijnw  o,  0*1,  130 
Nearr  \\  Point  04 
Far  the  :  t  Point  04 
Oja  I  --  see  pelnnnay  Diagram 
Spherical  73 

Edge  (Convey  Polygon)  76,  130 
SphcfK.al  130 
E clge  (Convc*  Polygon)  132 

Woide  11,  17,  107.  123,  133 

Yaglcui  and  BoltyaitAii  24 

Yao  0,  24.  101 

Yuval  3.  14,  17,  101.  113.  122 


Zolnowcky  7.  32.  127 
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